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Summary 


This report summarizes research work on the implementation of a partitioned numerical 
algorithm for determining the dynamic response of coupled structure/controller/estimator 
finite-dimensional systems. The partitioned approach leads to a set of coupled first and 
second-order linear differential equations which are numerically integrated with extrapola- . 
tion and implicit step methods. The present software implementation, ACSIS, utilizes par- 
allel processing te chni ques at various levels to optimize performance on a shared-memory 
concurrent/vector processing system. The current work also generalizes the form of state 
estimation, whereby the Kalman filtering method is recast in a second-order differential 
equation equivalent to and possessing the same computational advantages of the structural 
equations. As part of the present implementation effort, a general procedure for the design 
of controller and filter gains is also implemented, which utilizes the vibration character- 
istics of the structure to be solved. Example problems are presented which demonstrate 
the versatility of the code and computational efficiency of the parallel methods is exam- 
ined through runtime results for these problems. A user’s guide to the ACSIS program, 
including descriptions of input formats for the structural finite element model data and 
control system definition, can be found in Appendix A. The procedures and algorithm 
scripts related to gain design using PRO-Matlab axe included in Appendix B. In Appendix 
C, a stability analysis for the partitioned algorithm is presented which extends previous 
analysis to include observer dynamics, leading to a clearly definable stability limit. The 
source code for the parallel implementation of ACSIS is listed in Appendix D. 


1.0 Introduction 

The present work on the implementation of a partitioned transient analysis algorithm for 
the simulation of linear Control-Structure Interaction (CSI) problems has concentrated 
on four major areas. The initial software implementation emphasized the user-friendly 
aspect and a structural dynamics-oriented interface for experienced practitioners of finite 
element analysis progr ams . Another reason for a new implementation of the algorithm 
was that the initial architecture of the CS3 software testbed developed by Belvin and Park 
[1-2] had little provision for effective parallelization. CS3 also included extensive links to 
optimization and optimal control algorithms which were not central to the current work 
and proved to be a further hinderance. This new software implementation was designated 
ACSIS, for Accelerated Control-Structure Interaction Simulation, and led to significant 
improvements in speed for particular problems on conventional serial processors due to 
simpler and more economical storage of primary variables. ACSIS is also versatile in its 
usage of a general Timoshenko beam element with pin release capability and shear correc- 
tion factor adjustment, and control system definition via a single data file. A preliminary 
Users Guide was developed for ACSIS, and the input formats were made compatible with 
pre/postprocessing software developed for Sun and Silicon Graphics computers so that 


future versions using parallel techniques via element mesh domain decompositions can rely 
on the same X- Window-based I/O utilities. Table 1 documents runtime comparisions of 
ACSIS and CS3 on a sample 48 DOF problem simulated on a Sun 3/260 workstation using 
a floating point accelerator. 


2.0 State Estimation via Second-Order Kalman Filtering 

One restriction of the CS3 software testbed was the form of dynamic observer equations 
used in the partitioned algorithm as developed in [2]. However, Belvin and Park [3] showed 
that a general Kalman filtering type of state estimation was not only possible in the parti- 
tioned solution, but could be implemented at a very small additional cost in computations 
by a slight modification in the way the dynamical equations of the plant are cast into 
first-order form for filter gain design. The methods employed in [3] have been successfully 
implemented into ACSIS, thereby enhancing the code’s ability to handle a wide variety 
of state estimation schemes. This is important as the application of optimal control tech- 
niques to the state estimator problem leads to this more general form, and as the restric- 
tions of the observer used in CS3 typically meant either discarding part of the observer 
gain parameters, resulting in a loss in system performance. 

ACSIS retains the option of using the more resricted observer form, as well as simple 
full state feedback (no dynamic compensation). It is clear from examining the respective 
equations, that for structures with stiffness- proportioned damping, the only additional 
calculations required for the K alman filter is the multiplication of the Li gain matrix by the 
predicted state estimation error vector j (see equations 25 of [3]). This is not a significant 
portion of the computations required at each integration step, as the dimension of 7 is small 
(number of sensors). By far the major costs are for computing an internal force of the 
form Kq, and backsolving the factored integration matrix for the estimated displacement 
states. This is verified numerically in the ACSIS results of Table 1 (using a restricted form 
of observer) and version (Al) of Table 2. The model used in Table 1 is roughly comparable 
to the 54 dof truss in Table 2, and both show the additional simulation time due to state 
estimation is roughly the same as an additional transient analysis. Finally, the Kalman 
filtering equations do not lead to any additional complications in parallel implementation of 
the overall algorithm as compared to the more restricted second-order observer developed 
for CS3. 
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3.0 Parallel Implementation of ACSIS 

A primary emphasis in our work dealt with the optimization of the software implemen- 
tation on a concurrent processing system. The platform chosen (primarily due to avail- 
ability, initially at CU and later at NASA LaRC) was the Alii ant FX/8 shared memory 
multiprocessor system with 8 parallel processing units and vectorization capabilities. Ver- 
sions of the software have been ported to the Alliant and compiled using the Concentrix 
FX/Fortran optimizing compiler, which has available options for automatic vectorization 
and concurrency of standard, problem-independent parallel computations. 

The partitioned CSI algorithm has three primary levels of parallelism in its numerics which 
can be exploited. At the highest level is the integration of the second-order dynamical plant 
(structure) and filter (estimator) equations, which as designed are of roughly equivalent 
size. Through the algorithm, these systems are effectively decoupled and independent 
at each discrete time step, and thus may be handled in parallel by invoking a compiler 
directive in the main program, which calls the respective subroutines simultaneously and 
handles re-synchronization of the execution upon their return. 

At a lower level of the algorithm, the structure and state estimator equations exhibit 
symmetric, second-order forms typical of linear structured dynamics problems. It is well 
known that computations related to the formation of the internal force vector, Kq, can be 
re-implemented at an element level [4], which, through decomposition of the element mesh 
[5], can be handled in parallel within exclusive subdomains. This technique becomes par- 
ticularly attractive for larger problems on more massively parallel systems; in the present 
work, the element-by-element (EBE) technique was not as effective as other types of paral- 
lelism. It should be noted here that the partitioned algorithm employs implicit integration 
methods, leading to systems of algebraic equations which are factored and solved using di- 
rect, rather than iterative, numerical methods. Therefore, formulation of the internal force 
vector is needed only in the formation of the known (right hand side) vector for integration 
of the plant and filter equations. The alternative to EBE computations is multiplication of 
the relevant displacement vector with the global stiffness matrix stored in profile (skyline) 
form. To “simulate” the advantages of local memory typical in large-scale parallel process- 
ing systems, the computed element stiffnesses were saved in shared memory for the EBE 
calculations, thus avoiding the need to recompute this data at each integration step. In 
addition, a low-overhead automatic element domain decomposition was provided for the 
parallel EBE method. 

The final and lowest level of parallelism is obviously that of the basic matrix computations 
such as addition, multiplication, etc. These numerical operations me inherent in nearly 
all areas of the software implementation, including problem preprocessing. With the very 
capable optimizing compiler available on the Alliant system, this parallelism was exploited 
through vectorization and concurrency of the nominal source code using the FX/Fortran 
compiler run-time options -0 -DAS -alt. The performance of the resultant executable 
code was examined using the A lli ant’s profiling capabilities, and changes to the nominal 
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source code, remaining compliant to F77, to maximize the identifiable concurrency and 
thus enhance the resultant speed. This was particularly useful in the profile matrix/ vector 
multiply operation, whose speed is critical to the overall program performance. 


4.0 Problem Descriptions 


Three structural dynamics problems were developed for code testing at various levels of 
complexity. All three problems have the following common features: simulations consisted 
of 1000 integration steps and employed a stiffness-proportional damping. The damping 
was not needed for algorithm stability, but to ensure consistency between the examples 
and because the existence of damping in the plant equations has a strong influence on 
program speed. The Kalman filter models were of second-order form [2] and equivalent in 
size to their respective plant models. For controlled simulations, the control system began 
operating after 100 integration steps, and all gain matrices were full (i.e. all model states 
influence all actuators and are influenced by all sensors). Additional specific information 
for the problems follow. 


A. Axial Vibration of Elastic Bar (Spring Model) 


# Nodes: 

# Elements: 
# Degrees of freedom: 

# Actuators: 
# Sensors: 

Disturbance: 


3 free, 1 fixed 
3 
3 
1 
1 

Initial displacement 


B. Planar Vibration of Space Truss (Truss Model) 


# Nodes: 

# Elements: 
# Degrees of freedom: 

# Actuators: 
# Sensors: 

Disturbance: 


18 free 

33 

54 

4 

6 

Bang-bang type sinusoidal applied force 


C. General 3D Vibration of EPS Satellite Reboost (EPS7 Model) 


# Nodes: 

# Elements: 
# Degrees of freedom: 

# Actuators: 
# Sensors: 

Disturbance: 


97 free 
256 
582 
18 
18 

Bang-bang type square wave applied force 
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As can be seen, the problem sizes are roughly three different orders of magnitude, with 
corresponding increases in the sizes of the control systems. Appendix B includes routines 
using Matlab for the design of control and filter gains which were used for the control 
system design of all three example problems. An illustration of the EPS model is shown 
in Figure 1. 


5.0 Performance Assessment 

Table 2 compares CPU runtimes on the Alii ant computer (using the UNIX "time” com- 
mand) for distinct versions of the software. Version (Al) is the nominal F77 program code 
compiled without any performance-enhancing options, while version (A2) invokes auto- 
matic vectorization and concurrency of low-level, problem-independent computations such 
as vector addition and inner products. The performance improvements are significant, 
especially for the large EPS7 model, where the speed-up factor is 35-37. 

Version (A3) also uses the compiler options from (A2), but in addition has a compiler 
directive added to the main program which allows the plant and filter integration subrou- 
tines to be called in parallel. This does not affect the transient response results as that 
analysis option bypasses the altered code, but for controlled response there is some effect 
on performance. If filtering is used, which results in a significant increase in computation, 
the directive can lead to some increased speed as can be seen for the spring and truss prob- 
lems. There can also, however, be a reduction in performance as compared to (A2) if the 
finite amount of processors and vector units axe used in a less efficient way. This appears to 
be the case for the large EPS problem, where the "overhead” introduced by the directive, 
and its effect on processor assignment, is greater than the improvement generated by the 
manually-invoked parallel construct. 

Table 3 shows CPU run times for ACSIS using E-B-E computations, which, as mentioned 
previously did not lead to better performance on the example problems using the Alliant 
system. This appears to be due to the lack of effective vectorization of the individual 
element computations when forming the internal force via the EBE method. To deter- 
mine whether the EBE method effectively lead slower speeds through increased numbers 
of computations, version (A2) (see above) was altered by removing compiler optimization 
of the profile matrix/vector multiply routine (the alternate method to EBE). The resultant 
runtimes matched almost exactly with those of version (A4), leading to the conclusion that 
both methods require roughly the same amount of computations, but differ in how they 
can be optimized on the Alliant system. The matrix/vector multiply operation, in this 
environment, can exploit both vectorization and concurrency through the compiler’s per- 
formance options; this can be examined in the compiler output. The EBE computations, 
at the element level, do not vectorize because the parts of the global displacement and force 
vectors being operated on per element are not contiguous. In version (A5), the elements 
within each subdomain of the mesh are computed in parallel, leading to some performance 
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Overall, versions (A2) and (A3) provide the best code performance for the hardware avail- 
able. Parallelizing the observer and structure (A3) leads to mixed results; improvement 
for the small spring and truss problems, but not for the large EPS model. Element-by- 
element computations do not improve code performance over compiler optimization via 
vectorization and concurrency for this platform. Reimplementation of the algorithm lead 
to a 5:1 improvement over the CS3 testbed software on a serial computer (Table 1). Fur- 
ther optimization of ACSIS on the Alliant FX/8 lead to an additional 30:1 improvement in 
runtimes for large-order systems such as the 582 dof EPS model. Time history responses 
of selected variables for the example problems axe shown in Figures 2 through 10. For 
Figures 8 through 10, the u x ,u y ,u z displacements plotted axe located at node 45, which 
is located at the vertex of the large antenna of the EPS model (see Figure 1). 


6.0 Conclusions 

The present work has demonstrated the efficiency of a streamlined simulation code for 
the analysis of large-order CSI systems and the viability of the continuous second-order 
Kalman filtering equations for state estimation. These methods show the versatility of 
the partitioned CSI integration algorithm and the promise of its application to real-time 
simulation. It is evident, however, that the use of element-by-element computational tech- 
niques requires the development of innovative algorithms for effective implementation on 
massively parallel processing systems. Future work in this area will include integrating 
algorithms for on-line system identification and applying these capabilities to the problem 
of real-time control. 
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Simulation 

CS3 

ACSIS 

Transient 

439.2 

98.8 

Full State 

Feedback 

688.2 

181.5 

FSFB with 

Observer 

1156.7 

282.3 


Table 1: Comparison of Rimtime Speeds for CS3 and ACSIS 

on a Sun 3/260 System 




(Al) 

(A2) 

(A3) 


Problem 

Nominal 

Compiler 

Parallel 

Model 

Type 

Code 

Optimized 

Observer 

3 DOF 

Transient 

6.6 

2.1 

2.1 

Spring 

FSFB 

8.0 

3.3 

3.3 


K. Filter 

12.3 

3.5 

3.3 

54 DOF 

Transient 

78.2 

5.7 

5.6 

Truss 

FSFB 

97.1 

9.4 

10.2 


K.Filter 

170.7 

13.0 

10.7 

582 DOF 

Transient 

3506. 

98.6 

100.3 

EPS7 

FSFB 

7040. 

190.2 

294.5 


K. Filter 

n/a 

284.2 

312.5 


\ 


Table 2: CPU Results for Versions of ACSIS 
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Model 

Problem 

Type 

(A4) 

E-B-E 

Computation 

(A5) 

Parallel 

E-B-E 

(A6) 

Parallel 
Obs. & EBE 

3 DOF 

Transient 

3.8 

3.3 

3.3 

Spring 

FSFB 

4.9 

4.4 

4.9 


K. Filter 

6.6 

5.6 

5.0 

54 DOF 

Transient 

31.7 

13.0 

13.0 

Truss 

FSFB 

35.5 

16.9 

35.6 


K.Filter 

62.6 

27.3 

36.2 

582 DOF 

Transient 

391.7 

153.9 

n/a 

EPS7 

FSFB 

485.9 

245.9 

n/a 


K. Filter 

n/a 

n/a 

n/a 


Table 3: CPU Results for ACSIS with EBE Computations 
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Figure 1: EPS Finite Element Model 
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Truss Model: Open Loop Transient Response 
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Figure 5: Truss Transient Response 
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Truss Models Full State Feedback Response 
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Figure 6: Truss FSFB Response 
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Truss Model: Controlled Response w/Kalman Filter 
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Figure 7: Truss Response w/Filter 
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Accelerated 

Control 

Structure 

Interaction 

Simulation 


Introduction 

ACSIS is an analysis program for full-order simulation of control-structure interaction 
(CSI) problems. The CSI simulation is carried out using a partitioned analysis pro- 
cedure which treats the structure (or plant), the observer, and the controller/observer 
interaction terms as separate entities. This procedure allows ACSIS to maintain rela- 
tively small, sparse matrix equations when compared to the process of assembling the 
computational elements into a single set of equations of motion and solving simultane- 
ously. Although this software can also carry out modal analysis, the transient response 
and CSI simulation are done in real space using the entire finite element model. (For 
more information, see Ref. 2 of report, p. 14) 

The ACSIS program is run using two previously prepared data input files and interactive 
input of run options. The two data input files are the finite element input file and the 
controller definition file. 
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Interactive Options 

The run options are as follows, note that no defaults exist, data must be input each 
time it is requested for each item requested except during a background run. Files may- 
be given any names, example names are given only to match the truss example at the 
back of this manual. 

Please input analysis type: 

This option is for selecting modal analysis, CSI simulation or the transient response of 
a structure. Not ail interactive inputs are required for each analysis type. 

Do you wish to save an input file? (y or n) 

This option creates a file which saves all the interactive input options. ACSIS can be 
background run with minor option changes by editing this file and then directing the 
screen input to this file. This file is very different for each analysis type . (To do this 
run with the example names after running acsis.exe once interactively, the command 
would be: acsis.exe <INP_truss> &. 

Name of save input file? (filename) 

This question asks for the name under which to save the interactive input, example 
name: INP -truss 

Finite Element Model Input File Name (filename) : 

This file should contain all the finite element nodes, mesh, materials, properties, lumped 
inertias, fixations, and initial velocity and displacement conditions. It must be prepared 
in advance in the card format specified later, example name: FEM_truss 

Number of modes desired? - 

This only appears in the modal analysis to request the number of modes to be output. 
The actual analysis is carried out with double this number of modes for accuracy. 

Controller Definition File Name: 

This only appears for the CSI simulation. The file should contain the actuator and 
sensor locations, control gains, and observer gains. It also must be prepared in advance 
in card format, example name: C0N_truss 

Please input type of control:: 

This option is for selecting the form of control law equations used in the CSI simulation. 
Full state feedback uses the current states of the plant (structure) to determine the 
control via constant gain matrices. Second-order observer uses only the L2 filter gain 
matrix of a Kalman filtering type of state estimation where the state variables axe 
position and velocity. The Kalman filter option allows the use of a full set of filter 
gains, but the gain design must come from a alternate variable casting using position 
and generalized momenta. 
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Initial time, final time, control-on time, step size: 

Data format is four col umns for CSI, but only three columns for transient response since 
the control-on time part is deleted. 

Forcing function ID, scale factor, damping coeff- a,b: 

In order for any time dependent forcing function to be easily implemented despite vari- 
able time step sizes, the forcing functions must be entered into the subroutine f orces.f . 
The format is to assign each new forcing function an ID number in a elseif statement. 
An example of forces.f is included at the back of the manual. Then forces.f must 
be recompiled into an object file and acsis.exa relinked. This is easily done by typing 
make which detects changes to the fortran files and does only the necessary compiling. 
A particular forcing function is chosen by entering its ID number in the first column, 
in addition the force can be scaled by a constant using the scaling factor in the second 
col umn - The Rayleigh damping coefficients a and b are the third and fourth columns. 

Phase lag fix? (y or n) : 

Specifies whether to include an extra iteration of control and sensor state prediction 
at each integration step to improve accuracy. Not generally needed unless the user is 
investigating the source of instability in a simulation and needs to test the sensitivity 
of the response to the partitioned algorithm’s extrapolation method. 

Gain scale factors (4 toted) : 

These scaling factors are, in order: the FI control gain matrix (displacement), F2 control 
gain matrix (velocity), LI and L2 state estimator filter gain matrix. This question 
appears only for CSI. 

ACSIS has two types of output options. The first is the displacement or velocity motion 
of up to twenty separate degrees of freedom. Interactive questions are, ‘Number of 
displacement results to output (max 10)’ followed by as many ‘Input node #, 
dof for displacement output #’ as necessary. Then these repeat for velocity results. 
The name of the file where the output will be stored is requested by ‘Output file 
name? (filename).’ example name: OUT.truss. The second output option saves 
the displacement of all nodes at any time step where output is sent (see next entry). 
The format is suitable for animation of the entire structure. Question asks ‘Animation 
Output? (y or n) ’ then for an animation filename if necessary. 

Send output every how many steps? 

This option affects both output options to reduce the size of the output and animation 
files. It causes output to only be sent after a integer number of time step iterations. To 
get output each time step, simply enter 1. 
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Finite Element Input File 

The finite element input file consists o£ title cards followed by columns of data- The title 
cards can be in any order, but they must be all capitalized with the appropriate number 
of columns of data for each card. The program reads rows of data until encountering 
a new card. Data which represents an integer value may be entered with a decimal 
point while real data may be entered without a decimal point as necessary. Any line 
beginning with a * anywhere in the file is ignored and can be used to insert comments. 
Any blank line will result in a read error. 

NODES 

Each node must be defined on a separate line. Columns can be separated by any number 
of spaces or a comma. Data format is four.columns: node number, x-coordinate, y- 
coordinate, z-coordinate. 

TOPOLOGY 

Each element must be defined on a separate line. Truss elements require two nodes then 
two columns of zeroes. (Truss elements would also require pin releases in ATTRIBUTES 
below.) Beam elements require two nodes then a third reference node representing a 
point in the xz plane of the beam and then a column of zeroes. Element type refers to 
finite elment formalation. (Currently only type 1 = timeshenko beam element is now, 
available.) Data format is six columns: element number, element type, node #1, node 
■#2, node #3, node #4. 

ATTRIBUTES 

Each element is characterized by an ID number from each of the MATERIAL and PROP- 
ERTIES cards below. Each element also has six pin release codes. The first code is for 
longitudinal stiffness, the second code is for torsional stiffness, the third and fifth codes 
are bending stiffness at each end in the y direction and must be the same value, and the 
fourth and sixth codes are for bending stiffness in the z direction and also must be the 
same. (0 is stiff, 1 is released.) Data format is nine columns: element number, material 
type, property type, and six pin release codes. 

MATERIAL 

The material data is formatted in four columns: material type, Young’s modulus, shear 
modulus and density of the material. 

PROPERTIES 

The properties data is formatted in six columns: property type, cross sectional area of 
the element, I y * I z , I y , I z , shear shape factor SSF2, shear shape factor SSF3. 

FIXITY 

Nodes with any fixations are defined here in the finite element file. The nodes must be 
entered with the fixity of all six of their DOF’s, restrained or not. (1 is restrained, 0 is 
free) Data format is seven columns: node number, x, y, z, 4> z , <p y , <j> z . 
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INERTIA 


All lumped inertias must be entered with each separate DOF on an individual line. 
Therefore a single node could take up to six lines to define. Data format is three 
columns: node number, DOF number (1-6), and value of inertia. 

INITIAL CONDITIONS 

All initial displacement and velocity conditions are entered into the finite element file. 
Data format is four col umns : node number, DOF number (1-6), initial displacement, 
and initial velocity. 

END 

End of file. 

Controller Definition Cards 

The controller defintion file also consists of title cards followed by data entry. Each card 
must be followed by the appropriate number of columns of data and in some cases the 
appropriate number of rows. Integers can be entered in real format and vice versa if 
necessary. Any blank line will result in a read error. 

NACT 

Number of actuators in the entire control system. 

BMAT 

This entry creates the actuator position matrix or B matrix. There should be one 
row for each actuator, data format is four columns: node number, DOF number (1-6), 
actuator number, and sensitivity. 

NSEN 

i 

Number of sensors in the entire controls system. 

HDMA 

This entry creates the matrix of displacement sensor locations. One row per displace- 
ment sensor, data format is four columns: node number, DOF number (1-6), sensor 
number, and sensitivity. 

HVMA 

This entry creates the matrix of velocity sensor locations. One row per velocity sensor, 
data format is four col umns : node number, DOF number (1-6), sensor number, and 
sensitivity. 

F1GA 

This is a list of the Fl or displacement control gains. The data format is four columns: 
node number, DOF number (1-6), actuator number, and value of gain. 
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F2GA 


This is a list of the F2 or velocity control gains. The data format is four columns: node 
number, DOF number (1-6), actuator number, and value of gain. 

L1GA , 

This is a list of the state estimator Ll filter gains. The data format is four columns: 
node number, DOF number (1-6), actuator number, and value of gain. 

L2GA 

This is a list of the state estimator L2 filter gains. The data format is four columns: 
node number, DOF number (1-6), actuator number, and value of gain. 

END 

End of file. 

Examples 

This section includes all the files and procedures needed to run all three analysis on a 
simple elastic bar problem. The naming of the files is a simple and easy to remember 
system, however no particular format is necessary. 

The finite element file was created simply by typing the node locations, connectivity 
(topology), etc. with a text editor. 

File: FEM_spring 


MESH 

* 

NODES 

1 0.00 0.00 0.00 

2 1.00 0.00 0.00 

3 2.00 0.00 0.00 

4 3.00 0.00 0.00 

TOPOLOGY 

1 112 0 0 

2 1 2 3 0 0 

3 1 3 4 0 0 

ATTRIBUTES 

1 11011111 
2 11011111 
3 1 1011111 

MATERIAL 

1 1000 . 0.0 0.0 

PROPERTIES 
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1 1.00 0.00 0.00 0.00 0.0 0.0 

FIXITY 

1 111111 
2 0 11111 

3 0 11111 

4 0 11111 
INERTIA 

2 1 0.100 

3 1 0.100 

4 1 0.100 

INITIAL 

3 1 0.100 0.000 

END 


The modal analysis only requires the number of modes desired in addition to the finite 
element file. The file INP-springO documents the interactive inputs used in the modal 
analysis. The results axe saved in file EIG .spring. 


File: INP_spring0 


-1 

n 

* ACSIS input file, two lines above are 

* analysis type and save input file. Do 

* not change then by editing this file. 

* Finite element input file?(f ilenane) 
FEM_spring 

* Number of inodes desired? 

3 

* Output file? (filename) 

EIG.spring 


File: EIG_spring 


SUBSPACE ITERATION ROUTINE 

NB OF EIGENVALUES =» 3 

NB OF VECTOR* 3 

NB OF D0F» 3 
TOLERANCE- 1.000E-04 

NB OF RIGID MODES- 0 


ITERATION NO 1 
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l.OOOE+OO 1 . OOOE+OO l.OOOE+OO 

ITERATION NO 2 
1 .297E-14 3 . 194E-14 3.899E-14 

EIGEN ANALYSIS RESULTS: t < 

RADIAL - CYCLIC 

MODE EIGENVALUE FREQUENCY FREQUENCY 


1 

1980.6 

44.504 

7.0831 



2 

15550. 

124.70 

19.846 



3 

32470. 

180.19 

28.679 



EIGENVECTORS : 





1 

2.3305 

1.8689 

-1.0372 

0. 

0 . 

2 

1.8689 

-1.0372 

2.3305 

0 . 

0 . 

3 

1.0372 

-2.3305 

-1.8689 

0 . 

0 . 


MASS MATRIX DIAGONAL: 

213 l.OOOOOOOOOOOOOD-Ol 
312 l.OOOOOOOOOOOOOD-Ol 
411 l.OOOOOOOOOOOOOD-Ol 


To run the transient response of the structure, a forcing function or initial condition 
would be needed to excite the structure. An initial condition would be added to the 
finite element file. A forcing function must be added to f orces.f with a new ID number, 
then this number given as interactive input. In this case, an initial displacement acts 
on the second degree of freedom, which is defined in the finite element model input file. 
The file INP jspringl documents the interactive inputs used in the transient analysis. 

File: INP_springl 


-3 

n 

* ACSIS input file, two lines above are 

* analysis type and save input file. Do 

* not change them by editing this file. 

* Finite element input file? (filename) 
FEM. spring 

* Initial, final, step size? 
0.00000000 1.00000000 0.00100000 

* Forcing function, scale f, damping a,b? 
0 0.000000 0.00000000 0.00002000 

* Output file name? (filename) 

0UT_spring 

* Number of displacement outputs? 

1 

3 1 

* Number of velocity outputs? 
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0 

* Send output every how many steps? 

1 

* Send animation output?(y or n) 
n 

In addition to the finite element file, the interactive input, and an excitation, CSI 
simulation requires a controller definition file. A full state feedback controller for the 
truss structure is defined in COS_spring. 


File: COH_spring 


NACT 

1 

BMAT 

3 1 1 1.0 
F1GAIN 

211 193.31415000000 

311 1574.6315000000 

411 -1669.0460000000 

F2GAIN 

211 20.774256000000 

311 27.790403000000 

411 10.780212000000 

NSEN 
1 

HVMAT 
3 1 1 1.0 
L1GAIN 

211 8 . 7928909000000D-02 

311 -4. 880790 1000000D-02 

411 3 . 0102735000000D-03 

L2GAIN 

211 1 . 0380932000000D-0 1 

311 6.1410130000000 

411 -3.0608725000000 

END 



Then if acsis.exe is run interactively and the input is saved, the file INP_spring2 can 
be produced. The file could be edited to change the input files, the output file, the 
forcing function ED, the length of simulation, control-on time, etc. Then to run it again, 
type acsis.exe<INP_spring2> scr A. The scr is a scratch file which will store the 
screen output. 


File: INP_spring2 


2 
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* ACSIS input file, two lines above are 

* analysis type and save input file. Do 

* not change then by editing this file. 

* Finite element input file? (filename) 

FEM_ spring 

* Controller file name? (filename) 

CON. spring 

* Please input type of control: 

3 

* Initial, final, control-on, step size? 

0.00000000 1.00000000 0.10000000 0.00100000 

* Forcing function, scale f, damping a,b? 

0 0.000000 0.00000000 0.00002000 

* Phase lag fix?(y or n) 

n , 

* Gain scale factors (4 total)? 

1.00000000 1.00000000 1.00000000 1.00000000 

* Output file name? (filename) 

0UT_ spring 

* Number of displacement outputs? 

1 

3 1 

* Number of velocity outputs? 

0 

* Send output every how many steps? 

1 

* Send animation output? (y or n) 
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appendix b 


Matlab Procedure and Script, 
or ControUer and Kalman Filter 
/** * 
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Introduction 


In order to retain simplicity in the ACSXS code for parallel implementation purposes, 
no control system design algorithms were included. Instead, using modal data out- 
put from the eigenmode analysis module of ACSIS, a procedure was developed using 
the Pra-Matlab and its Control System Toolbox, which includes algorithm scripts for 
optimal control solutions via the solution of an algebraic Ricatti equation. In order to 
accomodate large order dyn ami cal systems, the design is accomplished in the uncoupled 
normal modes domain, using the available lowest eigenmodes from ACSIS. 

The procedure begins by copying and editing the mode data output from ACSIS (see 
the listing for EIG_spring in Appendix A) into readable variables for input to Matlab. 
The typical approach was to create one file as a Matlab script (i.e. a ” .m” file), with 
the eigenvalues, eigenvectors, and mass/dof data from ACSIS at the beginning, followed 
by actuator and sensor influence matrices (related to the physical degrees of freedom), 
the objective function weighting matrices, and the function which calls other Matlab 
scripts to determine the solution. An example of the above input for the spring problem 
described in Appendix A is in file EIG_spring.m, which is listed below. Compare this 
the the ACSIS mode data output shown in Appendix A to see how the editing was 
accomplished, and the additional control design data added 


File: EIG_spring . m 


lam ■ [ 

1980.6 

15550. 

32470. 

]; 




vi®[l 

2.3305 

1.8689 -1.0372 

0. 

0. 


2 

1.8689 

-1.0372' 2.3305 

0. 

0. 


3 

m*[ 2 

1.0372 
1 3 

-2.3305 -1.8689 

1 .00000000000003-01 

0. 

0. 

3; 


312 1.0000000000000D-01 

411 1.0000000000000D-01] ; 

t»[vl(:,2:4)3; 
qb«*zeros(3,l) ; 
qb (2,1) *1.0; 
hd=zeros(l ,3) ; 
h.v»zeros(l,3) ; 
hv(l,2)®1.0; 
qv»[l; .5; .1] ; 
q=diag( [qv ; qv] ) ; 
r*.0001*eye(l) ; 

[f l,f2]=alqr(lam,t,a,qb,q,r,3) ; 
qv*Cl;l;l] ; 
q*diag([qv;qv3) ; 
r*100*eye(l) ; 

Ckl,k23*akf (laa,t,a,hd,bv,q,r,3) ; 
save flout fl /ascii 
save f2out f2 /ascii 
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save klout kl /ascii 
save k2out k2 /ascii 


The scripts mlqr.m and mkf .m were written to accept as input the vector of eigenvalues, 
the eigenvector matrix (orthogonal vectors stored in columns), a matrix of node, com- 
ponent, d.o.f, mass data, and the actuator (or sensor) influence matrix and weighting 
matrices. The scripts output the gain results in four- column arrays, with one gain per 
row, and the corresonding node number, displacement component, and actuator (or 
sensor) identification. The top-level problem script (listed above) then saves the output 
in external files and the analysis is complete. The design scripts (listed below) also 
include checks on controllability and observability of the system based on the modal 
data and influence matrices defined, and produce plots of the closed loop poles resulting 
from the gain design. This aids the analyst in assessing the expected performance (and 
stability) of the exact system before moving the data back to ACSIS for simulation. 
The script contrank.m finds the rank of the controllability matrix through iterative 
rank calculations of submatrices so as to avoid the iilconditioning experienced in the 
full matrix. This is botfr.faster and more accurate for determining whether a particular 
actuator placement has full control of the included structural modes. 

File: mlqr.m 


f unct ion [? lout , F2out] =mlqr (lam , t , m , qb fl , R , nmode ) 

7 . 

7 . Controller gain design for second-order 
X structural system via given eigenmodes. 

X Gains are transformed to be coefficents 
X of structural variables (disp .velocity) ; 

X i.e. plant is second-order, of size ndof. 

7 . 

V, Arguments : 

7 . 

X lam: vector of eigenvalues (nmode z 1) 

X t: matrix of eigenvctors (ndof x nmode) 

X m: mass diagonal and dof mapping info (ndof x 4) 

X qb: actuator position influence matrix (nact x ndof) 

X Q: optimal design state weighting matrix (2*nmode x 2*nmode) 
X R: optimal design feedbk weight, matrix (nact x nact) 

X 

X Flout: FI gain matrix for 2nd-order plant 

X F2out: F2 gain matrix for 2nd-order plant 

7 . 

X Written by K.F. Alvin 
X 

format short e 
nmodmax“length(lam) ; 

[ndof , nact] s s ize (qb ) ; 

X 

X Variables : 

X mass: mass matrix 
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% A: state transition matrix 

% B: actuator influence matrix 

% G: control gain matrix 

X 

massd(m(: ,3))»m(: ,4) ; 
mass a diag(massd) ; 

A=» [zeros (nmode) .eye(nmode) ;-l*diag(lam(l :nmode)) .zeros (nmode)] ; 

Amax® [zeros (nmodmax) .eye(amodmax) ;-l*diag(lam) .zeros (nmodmax)] ; 

B* [zeros (nmode ,nact) ; t ( : , 1 : anode) 1 *qb] ; 

Bmax» [zeros (naodmax ,nact ) ; t ’ *qb] ; 

disp( ’Number of structural modes and actuators used: 1 ) 
disp ( [nmode ,nact] ) 

dispC’Rank of the controllability matrix: *) 
dispCcontrank(A.B)) 

disp (’Determining controller gains for given system. . . *) 

G«lqr(A,B.q.R); 

vmax a l.l*max(sqrt(laa)); 

axis ( [-smax , smax , -umax , smax] ) ; 

plot(eig(A-B a G) , ’♦*) 

grid 

title( ’Roots of controlled system*) 

hold 

pause 

x 

X partition gain matrix 

X 

G1 S G( : , i : nmode) ; 

G2=G( : , nmode* 1 :nmode+nmode) ; 

y. 

X Transform resultant gain matrices for use in 
X partitioned csi algorithm using second-order 
X structure (plant) equations. 

X 

disp (’Mapping ga:'ns bach to physical domain...’) 
fl a Gl*t(: ,l:nmode) ’♦mass; 
f 2=G2*t ( : , 1 : nmode) ’ ♦mass ; 

X 

X find modal damping ratios of controller for calulated gains: 

X 

Gmax»[f l^t ,f2*t] ; 
lambda a eig(Amax~Bmax*Gmax) ; 
plot (lambda , ’ ♦ ’ ) 
pause 

nfreq a sqrt(imag(lambda) ."2 + reed, (lambda) . *2) ; 
mdamp<*-real (lambda) ./nfreq; 

disp (’Resultant modal damping ratios for controller:’) 
disp([’ Damping *,* Damped Freq (rad/s) ’]) 
disp( [mdamp. nfreq. ♦sqrt (1-mdamp. ~2) .lambda]) 
bdamp a max(-2^mdamp . /nfreq) ; 

disp( ’Estimated mi ninn™ stiffness damping coefficient necessary’) 
dispC to stabilize residual modes due to gain roundoff accumulation: *) 
disp(bdamp) 

disp( ’Writing gains in node correspondence output form...’) 
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Flout a zeros(nact*ndof ,4) ; 

F2out=zeros (nact*ndof , 4 ) ; 
for i«l:nact; 
kmin 3 (i-l)*ndof+l; 
taaax 3 i*ndof ; 

Flout (taain:tanax, :) a [m(: ,1:2) ,i*ones(ndof ,1) ,f l(i,m(: ,3)) ’] ; 
F2out(taain:taaax, :) a [m(: ,1:2) ,i*ones(ndof ,1) ,f2(i,m(: ,3)) '] ; 
end; 

disp ( * Finished mlqr * ) 


File: mkf .m 


function IX lout ,L2out]*mkf (lam, t,m,hd,hv,Q,R, nmode) 

X 

X Kalman filter design for second-order 
X structural system via given eigenmodes 
X and transformed to independent 
X displacement/gen . momentum variable 
X casting for partitioned csi transient 
X analysis. See Belvin/Park paper for 
X filter variable definitions. 

X 

X Arguments : 

X 

X lam: vector of eigenvalues (nmode s 1) 

X t: matris of eigenvctors (ndof x nmode) 

X m: mass diagonal and. dof mapping info (ndof x 4) 

X hd: sensor position influence matrix (nsen x ndof) 

X hv: velocity position influence matrix (nsen x ndof) 

X Q: optimal design state weighting matrix (2*nmode x 2*nmode) 
X R: optimal design feedbk weight, matrix (nsen x nsen) 

X 

X Llout: LI gain matrix for 2nd-order filter 

X L2out: L2 gain matrix for 2nd-order filter 

X 

X Written by K.F. Alvin 

X 

format short e 
nmodmax a length(lam) ; 

[nsen, ndof] =size(hd) ; 

X 

X Variables : 

X mass: mass matrix 

X A: state transition matrix 

X G: noise influence matrix 

X C: output influence matrix 

X K: filter gain matrix > 

X 

massd(m(: ,3))=m(: ,4) ; 
mass=diag(massd) ; 

A 3 [zeros (nmode) ,eye(nmode) ;-l*diag(lam(l:nmode)) , zeros (nmode)] ; 
Amax 3 [zeros (nmodmax) , eye(nmodmax) ;-l*diag(lam) .zeros (nmodmax)] ; 
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G=eye(nmode+nmode) ; 

C*Chd*t(: ,l:nmode) ,hv*t(: »l:nmode)3 ; 

Cmax a [hd*t ,hv*t] ; 

disp (’Humber of structural modes and sensors used:*) 
disp ( Cnmode ,nsen] ) 

disp ('Rank of the observability matrix: ’) 
disp (eontr ank (A * , C * ) ) 

disp (’Determining filter gains for given system...*) 
K a lqeU,G,C,Q.R); 

X 

X partition gain matrix 
X 

Kl a K(l:nmode, :) ; 

K2=K (nmode+ 1 : nmode+nmode , ; 

Kmax a [K 1 ; zeros (nmodmax-nmode ,nsen) ; K2;zero3 (nmodmax-nmode , ns en)] ; 

X 

X find modal damping ratios of filter for calulated gains: 

X 

lambda a eig(Amax-Kmax*Cmax) ; 

plot (lambda, ’+*) 

hold 

pause 

al 3 ~real (lambda) ; 
bl a imag (lambda) ; 

disp (’Result ant modal damping ratios for filter:’) 
disp(C' Damping ’,’ Freq (rad/s) ’]) 
disp([al./sqrt(al.~2+bl.“2) ,bl]) 

X 

X Transform resultant gain matrices for use in 
X partitioned csi algorithm using second-order 
X Kalman filter approach. 

X 

disp( ’Mapping gains back to. physical domain...’) 
ll a t(: ,l:nmode)*Ki; 

12=mass*t ( : , 1 :nmode) *K2 ; 

disp( ’Writing gains in node correspondence output form;..’) 
Llout a zeros(nsen*ndof ,4) ; 

L2out a zeros(nsen*ndof ,4) ; 
for i a l:nsen; 
kmin a (i-l)*ndof+l; 
kmax a i*ndof ; 

Llout (fcmin : kmax , :) a Cm(: ,1:2) , i*ones (ndof ,1) ,ll(m(: ,3) ,i)] ; 

L2out (kmin : kmax , :)*Cm(: ,1:2) ,i*ones(ndof ,1) ,12(m(: ,3) ,i)] ; 
end; 

disp (’Finished mkf ’) 


File: eontr ank. m 


function maxrankacontrankCa.b) 
maxrank=0 ; 

[instate, nact} a size(b) ; 
i=0; 


36 





»* 


mat°b ; 

nevr ank=*rani (mat) ; 
while nesrank > maxrank 
maxrank»newrank ; 
i a i+l; 

mat* [mat , (a"i) *b] ; 
nevrank«rank(mat) ; 
if newrank a *nstate 
maxrank=newr ank ; 
and 

end 


Unfortunately, the external files created from Matlab with the gain results are written 
completely in terms of real numbers, whereas the first three columns are actually to be 
read by ACSIS as integers (they are used as indices). A separate utility was written to 
convert the format of these files; the source code is listed below. On Unix systems, the 
user simply assigns the standard input to be the current data file created by Matlab, 
and gives another file name for the standard output. The code is basically just a filter to 
chahge the three columns of indices to integers. The output cm then be pasted directly 
into the control definition file used by ACSIS. 


File: convcont .f 


program convcont 

parameter (NMAX a 100000) 
real*8 f(NMAX),v(4) 

integer node(NMAX) ,dof (NMAX) , act(NMAX) 
n=0 

100 read (*,*,err a, 200,end a 200) (v(i) ,i a l,4) 
n«n+l 

node (n) =>int (v ( 1 ) ) 
dof(n) a int(v(2)) 
act(n) =int(v(3)) 
f(n) »v(4) 
goto 100 

200 do 300 i»l,n 

print *, node(i) ,dof (i) ,act(i) ,f (i)' 
300 continue 

end 
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APPENDIX C 


Stability Analysis of a 
CSI Partitioned Simulation 
Algorithm with State Estimator 



The equation of the open-loop plant without passive damping in modal second-order 
form is , 

q + u 2 q = u (1) 

The controller uses a second-order observer to estimate the plant state, along with a 
full-state feedback control gain design. 

u = - ( rju 2 p + (up) 

p + u 2 p = u + ( 7 ^ 

7 -z-p 
z = q 


where q and p are the plant and estimator states, respectively, u is the control force, 7 
is the state estimation error, 2 is the sensor output, and 77, (,( are gain coefficients for 
position and velocity feedback, and the estimator filter. 

The partitioned analysis procedure uses a stabilized form of the control law and esti- 
mation error determination to reduce inaccuracies associated with the extrapolation of 
variables in the controller force prediction. A first-order filtering is achieved by taking 
the time derivative of (2a, c), 

u = — rju 2 p — (up 
j = 2 -p 


( 3 ) 


and then embedding the equations of motion through substitution for p. This leads 
to the following two'coupled, .first-order differential equations for the prediction of the 
control force u and state error 7. 





TJU 4 

0 


} p 


(*) 


Time discretization of ( 4 ) using an implicit midpoint rule leads to the following coupled 
difference equation: 


1 +• S(u 

6 


6((u 

1 + 6 ( 




J S(u 3 — T]U 2 
{ Su 2 


K‘-{ c r}p- 


( 5 ) 


where 6 = half-step size = j. Solving this equation requires knowledge of the plant (to 
obtain sensor output) and observer states. These values are extrapolated as: 


Using these equations to obtain u n+ » and 7 n+ j allows the plant and observer equations 
to be solved independently. Midpoint time integration of (1) and (2b) leads to the 
following equations: 

(1 + *V) <Z n+ * = + q n + Sq n 

(1 + SW) J> n+ * = 6 2 u n+ * +p n + Sp n + 6 2 £ 7 n+ * 

i n+ * = J («” +} - 5") 

P" +i -f(p" +i -P") ( 7 ) 

q n+1 = 2 q n+ * ~ q n 

p n+l = 2 q n+ $ - q n 

4 n+1 = 2 q n+ $ - q n 

p n+x = 2q n+ ± - q n 


Computational stability of the modal form of the CSI partitioned equations of motion 
using the aforementioned time discretization can be assessed by seeking a nontrivial 
solution of 



such that 

|A|<1 


(9) 


for stability. Subtituting (8) into (5-7), we obtain 


Jx = 0 

where 


Xi 

II 

ST 

+ 

u n+ i 

p n+ j 

p n+ % 

p" ?"} 

X 2 

= [*? +i 

7 n +i 

q n +i 

q n +» 

?” ?”] 


Ju J12 
_ J 21 ^22 


( 10 ) 


( 11 ) 
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10 0 

(rju 2 — S(u> 3 ) (1 + 6 £uj) 0 


0 -1 

0 0 


Jn = 


J 22 = 


0 



- 6 2 (1 + 6 2 u 2 ) 

0- 

- 1 

-6 

0 



0 

1 

“7 




1 

i 

7 

0 



0 

-2 


A + 1 

0 

0 



0 

0 



- 

-2 

0 




‘0 

0 0 0 

0 

O' 







0 

6(& 0 0 

0 

0 





J 12 = 


0 

—6 2 Z 0 0 

0 

0 






0 

O 

O 

KKJ> 

1 

0 

0 







0 

0 0 0 

0 

0 







.0 

0 0 0 

0 

0. 






- 

0 

0 0 0 

0 

0* 






o 

o 

<0 

1 

0 

1 



■w 



0 

o 

o 

n 

1 

0 

0 



J 

21 = 


0 

0 0 0 

0 

0 






0 

0 0 0 

0 

0 





• 

0 

0 0 0 

0 

0. 



' 1 

0 


0 

0 


0 


-1 ' 

-1 

(1 + so 

0 

0 


0 


0 

0 

0 


(1 + S 2 U 2 ) 

0 


- 

1 

-<S 

0 

0 


1 

7 

1 


1 

6 


0 

0 

0 


-2 

0 

A + l 

0 

. 0 

' 0 


0 

_2 


0 


A + 1 . 


-s 

C u 

0 

A + 1 


( 12 ) 


(13) 


(14) 


(15) 


A nontrivial solution to (10) is found from 

det J = 0 (16) 

which leads to the characteristic equation 

((1 - <5£) (1 + <5 3 Cw 3 - 6 2 tju 2 ) + <5f (1 + 8 2 u 2 ))z i 

+ (1 — 6£) -J- 6£ (l + S 3 ^u 3 — S 2 t 7u/ 2 ))z 3 

((1 — 6£) ( S 2 u> 2 + 6 2 tju> 2 ) + 6 2 u 2 (l + 6 3 £uj 3 — 8 2 rju 2 ) 

+6£ (6£u> + 8 2 u> 2 + 8 4 u 4 ))z 2 
(S 2 u> 2 (8 (u>) -f* 5£ ( 6 2 uj 2 4- S 2 rju> 2 ))z 

+S 2 u> 2 ( S 2 uj 2 + S 2 rju 2 ) = 0 


where 
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A = 


1 + z 
1 -** 


|A |<1 


Re(z ) < 0 


( 18 ) 


Thus, a test of the polynomial equation for possible positive real roots by the Routh- 
Hurwitz criterion indicates that the partitioned. approach as applied to the modal equa- 
tions give a computationally stable solution for no velocity feedback C = 0 provided 


h < 


y/fjU 


( 19 ) 




(20) 
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APPENDIX D 


ACSIS Source Code 

tf 
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/ 


File: Makefile 


.SUFFIXES : .f 

FFLAGS « 

.f .0: 

fortran 

.0 

-c $( FFLAGS) $*.± 



OBJS 3 

aesis.o 

aesisont.o 

addstf .0 

\ 


bcaa 3 d . 0 

forces .0 

input . 0 

\ 


pa&vmul . 0 

prepf ea.o 

profile.o 

\ 


read.o 

solver. 0 

aephlag.o 

\ 


zerovect.o 

Xu.o 

prepcon . 0 

\ 


control. 0 

secorder.o 

measure. 0 

\ 


eigens . 0 

singeig.o 

animout . 0 

\ 


stif f rc . 0 
kf liter. 0 

eatifvm.o 

renum.o 

\ 

acsis.«x«: 

S(OBJS) 



1 

fortran 

-0 acsis.exe 

$(FFLAGS) $(OBJS) 



File: shared. inc 


t 



c 

c — 

c 

C shared. inc (ACSIS database) 

C 

C — 

c 

c 

C Argument definitions: 


C 

ad amp: 

Rayleigh damping coefficient alpha 

C 

b: 

actuator location matrix (packed storage) 

c 

brow : . 

row number of corresponding real value in b 

c 

bcol: 

column number of corresonding real value in b 

c 

bval: 

number of nonzero values in b 

c 

bdamp: 

Rayleigh damping coefficient beta 

c 

confile: 

controller input file name 

c 

contype : 

id for type of control 

c 

coxy z : 

array of the x,y, and z components of each node 

c 

delta: 

one-half time step 

c 

delaq: 

one-half time step squared 

c 

ec: 

control prediction integration coefficient matrix 

c 

•mat: 

array of element material types 

c 

eo: 

observer construct matrix S in vector form 

c 

eprop : 

array of element property types 

c 

es: 

structure construct matrix S (H+delta«D+delsq*K) 

c 

etype: 

array of element types 

c 

elnum: 

array of element numbers for domain decomposition 

c 

f : 

vector of applied forces 

c 

fl: 

control gain matrix 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


f2 : 

t emf ile: 
forceid: 
gamma: 

S c: 

gk: 

go-' 

S» 5 

h: 

id: 

hdrow: 

hdcol: 

hdval: 

hv: 

hvroo: 

hvcol: 

hml: 

id: 

ix: 

inertia.: 

jdiag: 

11 : 

12 : 

mask: 

maaa: 

mat: 

mien: 

nact: 

ncsi: 

ndisont : 

ndof : 

ndomain : 

nal: 

neld: 

nap: 

naan: 

nvelout : 

nolag: 

outfile: 

out label: 

pin: 

pivot : 

prop: 

<15 

q.0: 

qalpha: 
qalphao : 
qbata: 
qbetao : 
qdot: 
qdotO : 
qe: 

qadot : 
r: 

scalef : 

stiff: 

tO: 

tc 

ti: 


control gain matrix 

finite element input file name 

identification number of forcing function 

atate correction force 

RHS vector for control prediction module 

SSS vector for Kalman filter momentum eqn 

RHS vector for obaerver module 

RHS vector for atructure module 

time atop size 

displacement sensor location matrix (packed storage) 

ros number of corresponding real value in hd 

column number of corresonding reed, value in hd 

number of nonzero values in hd 

velocity sensor location matrix (packed storage) 

rov number of corresponding real value in hv 

column number of corresonding real value in bv 

number of nonzero values in bv 

OOF slapping array: id(comp,node #)=Global DOF 9 

array of element connectivity and orientation 

array of concentrated inertias or lumped masses 

array of diagonal element addresses 

State estimator filter gain matrix 

State estimator filter gain matrix 

mass matrix N in reduced vector form 
array of different nmterials 

length of global matrices in profile vector storage 

actual number of actuators 

actual number of actuators and sensors 

number of displacement results to output 

actual number of degrees of freedom 

actual number of element domains (for dom. decomp.) 

actual number of elements 

array of number of elements in- each domain 

actual number of nodes 

actual number of sensors 

number of velocity results to output 

logical flag to signal corrector loop in measurement 

output file name 

array of output data requested 

array of element pin release codes 

Column pivoting info from FACTA 

array of different properties 

generalized displacement vector 

initial displacement condition 

gain scale factor for fi 

gain scale factor for 11 

gain scale factor for f2 

gain scale factor for 12 

velocity vector 

initial velocity condition 

state estimator displacement vector 

state estimator velocity vector 

Solution vector of control module {u > gamma}’ 

Scaling factor for forcing function 
stiffness matrix K in reduced vector form 
initial time 
control-on time 
final time 
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vector of control forces 


C 

C 

c 

c 

c 

c 

c 

c 

e 

e 

c 

c 


a: 


Parameter definitions 


MIX ACT: 
MAX CSX : 
MAID AT: 
HAXDOF: 
HAXBQM: 
HAXELE: 
HAXHLEH: 
HAXHODE: 


max. * of actuators 

max. combined 8 of actuators and sensors 
max. * of a&terials and properties 
aax. t of degrees of freed©® 
max. S of decomposition domains 
max. * of elements 

max. length of global vectors in reduced form 
sax. # of nodes 


parameter(HAXDOF=3000 , HAXACTsSO, HAXCSI=100) 
parameter(HAXH0DE=1000, HAXELE=3000, MAID 11=100) 
parameter (HAXHL£N=200000 , HAXD0H»50, HAXHZV=200) 

real* 8 to , tf , t c , h .delta , delsq , qalpha , qbeta , qalphao , qbetao 
real*8 q(HAXDOF) .qdot(HAXSOF) .qe(HAXDOF) .qedot(HAXDOF) 
real*8 n( MAX ACT) ,gamma(HAXACT) ,f (HAXDOF) ,r(KAXCSI) 
real*8 es(MAXMLEH) , eo(HAXHLEN) . ec(HAXCSI .HAXCSI) 
real*8 gs (HAXDOF) .go(MAXDOF) .gc(HAXCSI) .scalef ,pe (HAXDOF) 
real*? mass(HAXHLEH) ,stif (HAXHLEH) , ad amp, bdamp,gk( HAXDOF) 
r eal*8 eoxyz ( 3 , HAXHODE) .mat ( 6 , HAXD AT ) , prop ( 1 0 , MAXD AT) 
read *8 qO ( 6 , HAXHODE) , qdo 1 0 ( 6 , HAXHODE ) , inert ia ( 6 , HAXNQDE ) 
real«8 b(HAXHZV) .hd(HAXHZV) .hv(HAXHZV) ,estifm( 78.500) 
real*8 fi(HAXACT. HAXDOF) ,f2(HAXACT, HAXDOF) 
real*8 ll(HAXDOF.HAXACT) , 12 (HAXDOF .MAXACT) 
integer etype(HAXELE) , ix( 4, HAXELE) ,emat(MAXELE) .forceid 
integer eprop(HAXELE) ,pin(6 .HAXELE) ,id(6 .HAXHODE) 
integer mask(HAXHODE) .contype .brov(HAXHZV) .bcol(MAXHZV) 
integer hdrow(HAXHZV) .hdcol(HAXHZV) .hvrow(MAXHZV) 
integer hvcol(HAXHZV) ,bval,hdval,hvval 

integer ndof ,nact,nsen, ncsi, mien, jdiag(HAXDOF) .nnp.nel.neig 
integer neld(HAXDOH) .elnum(HAXELE.HAXDOH) , eldom(MAXELE) .ndomain 
integer outlabel(40) .ndisout .nvelout .pivot (MAXCSI) 
integer iadjcy (HAXHLEH) ,icount(HAXHODE+l) .perm(MAXHODE) 
integer xls (HAXHODE) 
logical animate, nolag 

character*32 f emfile.conf ile.outf ile.animf ile 


COHHOH 

CQHHOH 

COHHOH 

COHHOH 

COHHOH 

COHHOH 

COHHOH 

CQHHOH 

COHHOH 


COHHOH 


/FILES/ f emf ile , conf ile ,outf ile , outlabel , ndisout , 
nvelout , animf ile , animat e , no lag 
/TIHERS/ tO, tf.tc.h, delta, delsq 
/STATES/ q, qdot.qe.qedpt.u, gamma, f ,r,pe 
/FEHD AT/ mass , stif , adamp , bdamp , coxyz , mat , prop , qO , qdotO , 
inertia, scalef 

/IHTEGR/ es,eo,ec,gs , go ,gc,gk, pivot 

/DIKEHS/ ndof ,nact ,nsen ,nesi .mien , jdiag , nnp , nel ,neig 

/COHDAT/ b,hd,hv,fl,f2,ll,12. 


qalpha, qbeta, qalphao, qbetao 
/ELEDOH/ estifm,ndomain,neld,elnum,eldom 
/INTGEK/ f orceid , etype , ix , emat , epr op , pin , id , mask , 

contype , brow , bcol ,hdrow , hdcol , hvros , hvcol , 
b val , hd val , h wal 
/RESEQH/ iadjcy , i count , perm, r Is 


a 


File: acsis.f 
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CsPrograa ACS IS 

C=Purpoee Accelerated CSX Simulation 
C= Author K. Alvin 
C=Date Hay 1990 
CsBlock Fortran 


C — C 

C C 

C prograa ACSXS C 

C C 

C C 

C Purpose: 2nd Order Accelerated CSX Simulation C 

C C 

C C 

C — — C 

prograa ACSXS 

C GET SHARES DATA FILE 


include ’ shared . iac ’ 

C LOCAL VARIABLES 

real«8 t,z(MAXACT) 
integer n, a, runtype, outstrip 
C LOGIC 

call UiPUtCmntype , out skip) 

if (runtype) 100,200,300 

C EXGEHHODE AHALTSXS 

100 continue 

call PREPFEH 
call EIGEHS 
goto 999 

C CSX SIHULATIOH 

200 call PREPFEH 

call PREPCOH 

n a 0 
n a o 

print *, 'Finished Preprocessing . . . starting simulation’ 
print *,’Tiae »*,t0 
CALL ACSISOOT(tO) 

do 250 tatO.tf ,h 

call FQRCES(t*h/2) 

C Predict CSX coupling variables u and gamma 
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if (t . g« . tc) than 
call MEASURE (z) 
call COHTROl(s) 
if (aolag) then 
call SOPELA<I(*) 
call COHTROL(x) 
eadif 
•ad if 

C Structure and observer set ap for parallel execution 

CVS* CHCALL 

do 278 iel,2 

C Integrate Observer Equations 

if CCi .eq. 1) .and. (t .ge. tc)) then 
if (eontype . eq. 0) then 

Call SEC0RDER(aass , stif .adaap.bdamp.f ,go,eo ,qe,qedot 
delta , delsq , j diag , ndof , HAXDCF) 
elseif (eontype .eq. 1) then 
call KFILTER 
endif 

C Integrate Structure Equations 
elseif (i .eq. 2) then 

call SEC0RDER(aass , at if , adaap , bdaatp , f ,go , es , q , qdot , 
delta, delsq, jdiag, ndof , MAID OF) 

endif 

27S continue 

C PRIHT TIME EACH 100 iterations 

a * a + 1 ' 
a « a ♦ 1 

if (n .ge. 100) then 
print », ’Time ■ *,t+h 
n . » 0 
endif 

if (a .go. outskip) then 
call ACSISOUT(t*h) 

erite(24, ’ (40f 12 . 8) * ) t , (z(i) , i»l ,asen) 
a a o 
endif 

2S0 continue 
goto 999 

C TRAHSIEHT RESPOHSE 

300 call PREPFEM 

a a o 
a a o 

print «, ’Finished Preprocessing . . . starting simulation* 
print *,’Tiae a »,t0 
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call ACSISQtJT(tO) 
do 360 tetO,tf,h 
call F0ECES(t*h/2) 
call' ZEHQVECTCgs , ndof) 

call SECQEDE&Cmass , st if , adamp , bdamp , f , gs , e a , q , qdot , 
delta, delsq, jdiag.ndof .HAXDQF) 

C P&XST TIKE EACH 100 iterations 

a ■ a ♦ 1 
b ■ i ♦ 1 

if (a .go. 100) then 
priat e.'Time • ’ ,t*h 
a ■ 0 
eadif 

if (a .go. outskip) tkoa 
call ACSISOOT(teh) 
a ■ 0 
oadif 

3S0 coat inn o 

999 atop 
oad 


File: acsisout.f 


C=Hodule ACSIS0UT 

CaPurpose Vrito doairod output from ACSIS for plotting, otc. 

C=> Author K. Alvin 
C=Date Hay 1990 
C=Block Fortran 

C 

C Subroutine ACSISOUT 

C 

C Purpose: 

C This subroutine outputs formatted displacement and velocity 

C results at a given time for plotting time histories. The 

C desired output variables are defined in outlabel(). 

C 

C 

C Arguments 

C t - time 

C 


subroutine ACSISQUT(t) 

include ' shared. inc* 
real*8 t 
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c LOCAL VARIABLES 

iatagar i 
C LOGIC 

writa(i3, ’ (40212.8) ’ ) t,(q(id(ottlabal(i*10) ,oatlab«l(i))) , 

. i»i,adiaoot) ,(qdot(id(oatlabal(i+30) e oatlabal(i*20)) ) , 

. i»l,aoaloat) 

wrlta<a3,'(40 1*12.8)’) t- > (q a C i d(aat 3a b alCi»lft)«aatlabalCijjJ < 
. i»l „ndiaoat) , (qadot(id(oatlabel(i+30) ,oatlabal(i+20))) , 

. i*l,avalout) 

vrita(2S,’ (40*12.8)’) t,(a(i) ,i=l,nact) ,(gamma(i) ,ial,naea) 

i* (aaiaata) call AHIM0UT(q,id,aap,t,18) 

return 

«od 


File: addstf.f 


— +c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C 

C ARGUMENTS 

r«al*8 sk(naaq,asaq) , bk(l) 
iatagar ln(18) , jdiag(l) 

C iatagar 1 b( 18)7 jdiag(l), asaq 

C LOCAL ARGUMENTS 

iatagar i, j, k, 1, a 

C ASSEMBLE GLOBAL STIFFNESS AND LOAD ARRAYS 

do 20 j a 1, aaaq 
k « la(j) 

if (k . aq. 0) goto 20 


CsModala ADDSTF 

Caparposa Aaaambla Global ati**naaa matrix 
CaAatbor oho kaooa 
CaUpdata January 1989, by E. Praaoao 
C=Block Fortran 

subroutine ADDSTF(ak,ln,bk, jdiag.aaeq) 

C PURPOSE: 

C THIS SUBROUTINE ASSEMBLES THE ELEMENT STIFFNESS MATRICES 

C INTO THE COMPACTED GLOBAL STIFFNESS VECTOR. 

C 

C ARGUMENTS: 

C ak ELEMENT STIFFNESS MATRIX 

C la LOCATION VECTOR FOR ELEMENT STIFFNESS MATRIX 

C bk COMPACTED GLOBAL STIFFNESS VECTOR 

C jdiag - VECTOR OF DIAGONAL ELEMENT ADDRESSES 

C aaaq - NUMBER OF DEGREES OF FREEDOM PER ELEMENT 
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► 


1 * jdiag(k) - k 

C la jdlag(k*l) - k 

do 10 i » 1, aseq 
n a la(i) 

if(a .gt. k .0E. a .eq. 0) goto 10 

a ■ 1 ♦ a 

bk(a) a bk(a) ♦ ak(i,j) 

10 continue 

20 continue 
C 

return 

aad 

CaEnd Fortran 


File: beam3d.£ 


C=Module BEAM3D 

C=Purpose Conatract 3-d Tiaoahanko baaa aleaant stiffness aad lamped mass 
CaAuthor K. Alvin 
C=Data Nay 1990 
CaBlock Fortran 

sabroatina BEAH3D (n , ni , a j , ak , xyz , eaod , gaod , rho ,araa,ssf2,ssf3 t 
tc jtor,i2,i3, ipia.sk, am) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


ARGUMENTS : 

a Elaaant ID Number 

ni Hoda ID Number at End i 

nj Node ID Number at End j 

xyz Node Location Array 

amod Material Elastic Nodalaa (Young's Modulus) 

gmod Material Modulus of Rigidity (Shear Modulus) 

rho Material Mass Density 

area Element Cross-sectional area 

saf2 Shear shape factor in element x2 direction 

ssf3 Shear shape factor in element x3 direction 

jtor Torsional constant J 

12 Area moment of inertia about element x2 axis 

13 Area moment of inertia about element x3 axis 

ipin Pin release codes: 0»Fixed, l=Freed 

(1) Axial 

(2) Torsional 

(3) End A rotation about x2 axis 

(4) End A rotation about x3 axis 

(5) End B rotation about x2 axis 

(6) End B rotation about x3 axis 

sk Element Stiffness Matrix 

sm Element Mass Matrix 

integer n,ni,nj ,nk,ipin(l) 

real*8 xyz(3,l) ,emod t gmod,rho,area,ssf2,ssf3, jtor,i2,i3 
real*8 sk(12,l) ,sm(12,l) 

LOCAL VARIABLES: 

integer i.j 

real*8 dc(3,3) , length , rlength ,kc( 10) ,mc(3) 
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C LOGIC 


C Find El— at Length 

length * O.OdO 
do 10 i * 1,3 

de(i ,i) s xyz(i.nj) - xyz(i,ni) 

Xuagth s length ♦ dc(l,i)**2» 

10 continue 

-length » a qri (length) 
if (length .eq. O.OdO) then 

print *, ’3AH2D: Zero element length; n<* * ,n 

return 

endif 

C Find direction coeinee for xi,x2,x3 element axee 
do IS i*l,3 

dc(l,i) a de( 1,1) /length 
if (nk .eq. 0) then 
de(2,i) » 0.0 
else 

dc(2,i) * xyz(i.nk) - xyzCi.ai) 
endif 

15 continue 

if (nk .eq. 0) dc(2,3) » 1.0 

16 dc(3,l) ■ dc(1.2)*dc(2.3) - dc(2,2)*dc(i,3) 

dc(3,2) ■ dc(2,l)*dc(l,3) - dc(l,l)*dc(2,3) 

dc(3,3) » dc(l »l)*dc(2,2) - dc(2,l)*dc(l,2) 

rlength « eqrt(de(3,l)*»2 ♦ dc(3,2)**2 ♦ dc(3,3)**2) 

if (rlength .ne. 0.) goto 17 

dc(2,2) ■ 1.0 

dc(2,3) ■ 0.0 

goto 16 

17 do 18 i»l,3 

dc(3,i) « dc(3,i) /rlength 

18 continue 

dc(2,l) * dc(3,2)*dc(l,3) - d'c(l,2)*dc(3,3) 

dc(2,2) » dc(l,l)*dc(3,3) - dc(3,l)*dc(l,3) 

dc(2,3) « dc(3,lWc(l,2) - dc(l,l)*dc(3,2) 

C Compute various stiffness constants, accounting for pin codes 

if (ipin(l) . eq. 0) then 
kc(l) ■ are&*emod/length 

else 

kc(l) a O.OdO 
endif 

if (ipin(4) .eq. 0) then 
if (ipin(6) .eq. 0) then 

kc(2) a area*gmod*Bsf 2/length 
kc(6) a i3*emod/length 
ke(7) a kc(2)*length/2.0d0 
kc(9) a kc(7) *length/2 . OdO 
else 

print *,'BEAH3D: Pin code error, x3 direction, el #’,n 
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endii 


alaa 

it (ipia(6) . aq. 0) than 

print * , ’ BEAM3D : Pin coda arror, x3 direction, al #' ,n 

alsa 

kc(2) a O.OdO 
kc(6) * O.OdO 
kc(7) a O.OdO 
kc(9) • O.OdO 
andil 
+ndit 


it (ipin(3) .aq. 0) than 
ii (ipin(S) .aq. 0) than 
kc(3) * are**gBod*aai3/lehgth 
kc(5) a i2*eaod/length 
kc(8) a kc(3)*length/2.0d0 
kc(lO) a kc(8)«length/2.0d0 
alaa 

print *,’BEAH3D: Pin coda arror, x2 direction, al t ’ ,n 
eadii 

alaa 

ii (ipin(5) .aq. 0) than 

print *, ’BEAH3D: Pin coda error, x2 direction, al 8* ,n 

alaa 

kc(3) • O.OdO 
kc(5) a O.OdO 
kc(8) a O.OdO 
kc(lO) a O.OdO 
endif 
endif 

ii (ipin(2) .aq. 0) than 
kc(4) ■ jtor*gBod/length 

alaa 

kc(4) » O.OdO 
endi t 

mc(l) a area«rho*laagth/2 . OdO 
mc(2) a i2*rho*length/2.0d0 
mc(3) a i3*rho*length/2.0d0 


ak(l.l) 

sk(l,2) 

ak(l,3) 

sk(l,4) 

akCl.S) 

ak(l,6) 

ak(l,7) 

ak(l,8) 

ak(l,9) 

ak(l,10) 

ak(l.U) 

ak(l,12) 

sk(2,2) 


» kc(l)«dc(l,l)*dc(l ,1) 
kc(3)*dc(3,l)*dc(3,l) 
» kc(i)*dc(l,l)*dc(l,2) 
kc(3)*dc(3,l)*dc(3,2) 
a kc(l)*dc(l ,l)*dc(l,3) 
kc(3)*dc(3,l)*dc(3,3) 
» kc(7)*dc(2, l)*dc(3,l) 
a kc(7)*dcC2,i)*dc(3,2> 
a kc(7)*dc(2,l)*dc(3,3) 
a -ak(l,l) 
a -ak(l,2) 
a -ak(l,3) 

* ak(l,4) 
a ak(l.S) 
a ak(l,8) 

a kc(l)*dc(l,2)*dc(l,2) 
kc(3)*dc(3,2)*dc(3 ,2) 


♦ kc(2)*dc(2 1 ,l)*dc(2,i) * 

♦ kc(2)*dc(2,l)*dc(2,2) ♦ 

♦ kc(2)*dc(2,l)*dc(2,3) ♦ 

- kc(8Wc(3,l)*dc(2,l) 

- kc(8)*dc(3,i)*dc(2,2) 

- kc(8)«dc(3,l)*dc(2,3) 


♦ kc(2)*dc(2 ,2)»dc(2,2) ♦ 
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ak(2»3) ■ ke(l)*dc(l p2)*dc(l ,3) * kc(2)*dc(2,2)*dc(2.3) ♦ 

kc(3)*dc(3,2)*de(3,3) 

ak(2,4) » kc(7)*dc(2,2)*de(3,l) - kc(8)*dc(3,2)*dc(2,l) 

sk(2,5) « kc(7)*dc(2.2)*de(3,2) - kc(8)*de(3,2)*de(2,2) 

#k(2,8) » ke(7)*dc(2,2)adc(3,3) - kc(8)*de(3.2Wc(2,3) 

#k(2,7) a -sk(i,2) 
ak(2,8) a =sk(2,2) 
ak(2,9) a ~ak<2,3) 

«k(2,10) > ak(2,4) 

ak(2.1i) a ak(2,8) 

»k<3,12) « ak(2,6) 

ak(3,3) ■ kc(i)*dc(l,3)*dc(i,3) <• ke(2)*dc(2,3)*dc(2,3) ♦ 

ke(3)«de(3,3)«de(3,3) 

ak(3,4) a kc(7)*dc(2,3)*dc(3 ,1) - kc(8)*dc(3,3)*dc(2,l) 

sk(3 ,5) a ke(7)*de(2,3)*dc(3,2) - kc(8)*dc(3,3)*dc(2,2) 

ak(3,8) a kc<7)*dc(2,3)*dc(3,3) - kc(8)*dc(3,3)*dc(2,3) 

ak(3,7) » -ak(l ,3) 

sk(3,8) « -ak(2,3) 

ak(3,9) a -sk(3,3) 

sk(3p J.0) a ak(3,4) 

sk(3,ll) a sk(3,S) 

ak(3,12) a S k(3,6) 

8k(4,4) a kc(4)*dc(l ,l)*dc(l ,l)+’(kc(10)+kc(5)) , i'dc(2 ,l)*dc(2 , l) 

♦ (kc(9)*kc(6))*dc(3 , l)*dc(3,l) 

sk(4 ,5) a kc(4)ade(i,i)adca,2)+(kc(10)+kc(S))*dc(2,l)*dc(2,2) 

♦ (kc(9)+kc(6))*dc(3 ,l)*dc(3,2) 

ak(4,6) a kc(4)*dc(l,l)*dc(l,3)*(kc(10)*kc(5))*dc(2.1)*dc(2,3) 

♦ (kc(9)*kc(6))*dc(3,l)*dc(3,3) 

■k(4,7) a -«k(l,4) 

sk(4,8) » -sk(2,4) 

ak(4,9) a -ak(3,4) 

sk(4,10) a -ke(4)*dc(l,l)adc(l,l)*(kc(10)-kcC5))*dc(2,l)*dc(2,l) 

♦ (kc(9)=>kc(6))*dc(3,l)-*dc(3,l) 

sk(4,li) a -kc(4)*dc(l,l)»dc(l,2Wke(10)-kc(5))*dc(2,l)*dc(2,2) 

♦ (kc(9)-kc(8))*dc(3,l)*dc(3,2) 

ak(4,12) a -kc(4)adca,l)*dc(l,3)*(kc(10)-kc(5))*dc(2,l)*dc(2,3) 

+ (kc(9)-kc'(6))«dc(3 , l)*dc(3 ,3) 

sk(5,S) » kc(4)adc(l,2)*dc(l,2)+(ke(10)*kc(S))»dc(2,2)*dc(2,2) 

♦ (kc(9)*kc(6))*dc(3 ,2)*dc(3 ,2) 

sk(5,6) a kc(4)*dc(l,2)*dc(l,3)*(kc(10)*kc(5))*dc(2,2)*dc(2,3) 

♦ (kc(9)*kc(6))*dc(3,2)*dc(3,3) 

sk(5,7) a -ekU.S) 

ak(B,8) a -ak(2,S) 

sk(B ,9) a -ak(3,5) 

ak(S,10) a ak(4,ll) 

akCS.ll) a -kc(4)*dc(l,2)*dc(l ,2)*(kc(10)-kc(5))*dc(2,2)*dc(2,2) 

♦ (ke(9)-kc(6))*dc(3 ,2)*dc(3 ,2) 

ak(5,12) a -kc(4)*dc(l ,2)*dc(l ,3)+(kc(10)-kc(5))*dc(2 ,2)*dc(2,3) 

♦ (kc(9)-kc(6))*dc(2,3)*dc(3 ,3) 

ak(S,6) a kc(4)*dc(l,3)*dc(i,3)*(kc(10)*kc(5))*dc(2,3)*dc(2,3) 

♦ (kc(9)+kc(6))*dc(3,3)*dc(3,3) 

ak(6,7) a -ak(l,6) 

ak(6,8) a -ak(2,6) 

ak(6 ,9) a -sk(3,6) 

ak(6,10) a ak(4,12) 
ak(6,ll) a ak(5,12) 

ak(6,12) a -kc(4)*dc(l,3)*dc(l,3)+(kc(i0)-kc(5))*dc(2,3)*dc(2,3) 

♦ (kc(9)-kc(6))*dc(3 ,3)*dc(3 ,3) 
ak(7,7) a ak(l.l) 
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ak(7.8) a ak(l,2) 
ak(7.9) « ak(l,3) 

•k(7,10) - -ak(l,4) 

ak(7,U) ■ -ak(l,5) 
ak(7 , 12) ■ -ak(l,8) 

ak(8„8) ■ ak(2,2) 

ak(8,9) > ak(2,3) 

ak(8,10) > -ak(2,4) 

sk(8,lt) » -ak(2,6) 
ak(8,12) • -ek(2.6) 

ak(9,9) > ak(3,3) 

ak(9.10) » -ak(3,4) 

ak(9,U) a -ak(3,6) 
sk(9,12) > -ak(3 ,6) 

3k(10,10) « sk(4,4) 

ak(lO,ll) > ak(4,6) 
ak(10,12) - ak(4,6) 
akCll.ll) a ak(S.S) 
ak(ll,12) a ak(5,6) 
ak(l2,12) > »k(8,6) 

do 60 i a 1,12 
do SB j • 1,12 
an(i,j) > o.dO 
55 continue 

SO continue 



C Row-sum rotated mass matrix to re-diagoaalize 


C 

C 

C 

C 

C 

C 


C 

C 

C 

C 

C 

c 


sm(l,l) a mc(l) 

am(2,2) a mc(l) 

sm(3,3) a mc(l) 

am(4,4) a mc(2)*(dc(l,l)*dc(l, l)+dc(2,l)*dc(2,l)) ♦' 

»c(3)*(dc(l ,l)*dc(l,l)+dc(3,l)*dc(3 , 1)) ♦ 
mc(2)*(dc(l ,l)*dc(l ,2)*dc(2,l)*dc(2,2) ) ♦ 

. mc(3)a(dc(l;l)*dc(l,2)-*dc(3,l)*dc(3,2)) ♦ 

mc(2)*(dc(l ,l)*dc(l ,3)+dc(2,l)*dc(2 ,3) ) ♦ 

mc(3)*(dc(l ,l)*dc(l,3)+dc(3,l)*dc(3,3)) 
am(4,S) a mc(2)*(dc(l ,l)*dc(l ,2)+dc(2,l)*dc(2 ,2) ) ♦ 

ac(3)*(dc(l,l)*dc(l,2)+dc(3,l)*dc(3,2)) 
am(4,6) a oc(2)e(dc(l,l)*dc(l,3)+dc(2,l)*dc(2,3)) ♦ 

mc(3)*(dc(l,l)*dc(l ,3)+dc(3,l)*dc(3,3)) 
8a(5,4) a oc(2)*(dc(l,l)*dc(l,2)«>dc(2,l)*dc(2,2)) ♦ 

mc(3)*(dc(l,l)*dc(l.2)t'dc(3,l)*dc(3,2)) 
80(5,5) a oc(2)*(dc(l,2)*dc(1.2)+dc(2,2)*dc(2.2)) ♦ 

oc(3)*(dc(l,2)adc(l,2)*dc(3,2)»dc(3,2)) ♦ 

mc(2)*(dc(l,l)*dc(l ,2)*dc(2,l)*dc(2,2)) ♦ 

mc(3)*(dc(l,l)*dc(l,2)+dc(3,l)*dc(3,2)) ♦ 

oc(2)*(dc(l,2)«dc(l ,3)+dc(2,2)*dc(2,3)) ♦ , 

oc(3)*(dc(l ,2)*dc(l ,3)*dc(3 ,2)*dc(3 ,3)) 
ao(5,6) ■ ac(2)*(dc(l ,2)*dc(l,3)+dc(2,2)*dc(2,3)) ♦ 

oc(3)*(dc(l ,2)*dc(l ,3)*dc(3,2)*dc(3 ,3)) 
ao(6,4) » oc(2)*(dc(l,l)adc(l,3)+dc(2,l)*dc(2,3)) ♦ 

mc(3)*(dc(l , l)*dc(l ,3)+dc(3, l)*dc(3,3)) 
ao(6,5) ■ mc(2)*(dc(l ,2)*dc(l ,3)+dc(2,2)*dc(2,3)) ♦ 

oc(3)*(dc(l,2)*dc(l,3)+dc(3,2)*dc(3,3)) 
ao(8,6) a oc(2)*(dc(l,3)*dc(l,3)+dc(2,3)*dc(2,3)) ♦ 

oc(3)*(dc(l,3)*dc(l,3>dc(3,3)*dc(3,3)) ♦ 
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sa<7,7) 
sa(8,8) 
sa(9,9) 
sh(10,10) 
C sa(iO„ii) 

C s»(10,12) 

C sa(ll,l 0) 

«*(11 9 4D 
C s®(ii,12) 

C sb( 12,10) 

C sa(12,li) 

sa(12,12) 


ac(2)*<de(i,l)*dc(i,3)*dc(2,l)*de(2,3)) + 

ac(3)*<dc<l,l)*dcU,3)+dc(3,l)*dc(3 9 3)) ♦ 

■c(2)*(de(l,2)*dc(l,3)+dc(2,2)«dc(2,3)) * 

ac(3)*(dc(i,2)*de(l,3)+dc(3,2)*de(3,3)) 

ae(i) 

ac(i) 

ae(l) 

sa(4,4) 

s»(4„6) 

sa(4,6) 

ss(5,4) 

.sa£&*S) 

sa(f>,6) 

sa(6,4) 

sb(6,S) 

sm(8,6) 


do 100 i»i,12 
do 200 j»l,i-l 
sk(i, j) a sk(j,i) 
200 continue 

100 continue 


return 

end 

C=End Fortran 


File: forces. f 


CaKodule FORCES 

C=Purpose Calculate applied force rector at given time 
C=Autho r K. Alvin 
C=Date Hay 1990 
C=Block Fortran 

C 

C Subroutine FORCES 
C 

C Purpoae : 

C Returns force from stored function at any given time. 

C The forcing functions are hardwired by the user. The 

C function is selectable at program ezcecution using the 

C forcing function 10, which by convention is the statement 

C label used in branching. 

C 

C ; 

c 


subroutine FORCES (time) 

include • shared . inc ' 
real*8 time, pi 

C LOQXC 

pi > 3. 141S926 
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call ZEHOVECKf.ndof) 

101 if (forcaid . aq. 101) thaa 

if (tima .la. .02) thaa 

f (id(2,15)) ■ 100 . a ( 1 . -coa ( 2 . *pi*t ima/ . 02) ) 
a ad if 

102 alaaif (forcaid . aq. 102) thaa 

if (tima .It. .1) than 
f(id(2.96)) ■ 10. 
alaaif (tima .aq. .1) than 
f (id(2,95)) * 0. 

alaaif ((tima .gt. .1) .aad. (tima .It. .2)) thaa 
f(id(2,96)) « -10. 

alsa 

f(id(2,9S)) a o. 
aadif 

103 alaaif (forcaid .aq. 103) thaa 

if (tima .la. .01) thaa 

f (id(l.lS)) a ioo 
aadif 

104 alaaif (forcaid .aq.104) thaa 

if ((tima .gt. 0) .aad. (tima .It. .17)) thaa 
f(id(3,126)) a io 

alaaif ((tima .gt. .17) .aad. (tima .It. 1.0)) thaa 
f(id(3,12S)) a -io 
aadif 

105 alaaif (forcaid .aq. 105) than 

if (tima .la. .01) thaa 

f(id(2,9)) * I00.*(l.-coa(2.*pi*tima/.01)) 
alaaif (tima .la. .02) thaa 

f(id(2,9)) ■ 100 . *( cos (2 . *pi*tiaa/ . 01 )-l . ) 
aadif 

aadif 

do 10 i a i, adof 

f(i) * acalaf * f(i) 

10 coatinaa 

ratarn 

aad 


File: input. f 


C=Hodula IS POT 

C=Puxpoaa Input data paramatara for ACSIS 
C= Author K. Alvin 
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CsDate Hay 1990 
C=Block Fortran 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


Subroutine INPUT 


Argument definitions 


runtype 

savin 

runfile 

COBBSnt 

outskip 


10 of analysis run type 

variable to control creation of input file 
variable indicates if run is froa input file 
dummy naae for comment input lines 
nuaber of steps to skip before sending output 


subroutine INPUT (runtype, out skip) 


include ' shared. inc' 
integer runtype, outskip 
character 1 *! savin, ruaf ile, temp 
char act er *48 comment , inpf ile 


C PRINT AND READ START-UP 


print *, ’2nd Order Accelerated CSI Simulation (ACSIS)' 
print * 

print *, 'Please input analysis type:* 
print * 

print «, * 1. Eigenmode Analysis' 

print *, ’ 2. CSX Simulation' 

print *, ' 3. Transient Response' 

print • 

read ■*, runtype 


C RUN OPTIONS AND INPUT FILE SETUP 


runfile ■ 'n' 
if (runtype .It. 0) then 
runtype -1 * runtype 
runfile » 'y' 
endif 

print *, 'Do you oish to save an input file? (y or n)’ 
read 21, savin 

20 format (a32) 

21 format (ai) 

if (savin .eq. * y ’ ) then 

print *, ’Name of save input file? (filename)’ 
read 20, inpf ile 
open( 16 , f ileainpf ile ) 
runtype a -1 *• runtype 
vrite(16 , ’ (i2) 1 ) runtype 
. runtype » -1 * runtype 
vrite(16, ’(al)*) ’n' 

write(16,’(a47)’) '* ACSIS input file,tvo lines above are' 
vrite(16,’(a48)') '* analysis type and save input file. Do' 
vrite(16, ’ (a48) ’ ) ’* not change them by editing this file.’ 
endif 
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ruatype a runtype - 2 

if (runfile .aq. ’y’) than 
do 30 i a 1,4 

road 20, comment 
30 eontinno 

endif 

print *, ’Finite Element Xodal Input File Name (filename)’ 
road 20, femfile 
open (11 ,f ile=f oaf ilo) 
if (savin .aq. ’y’) than 

writo(18, ’ (a47) ’ ) ’* Finite element input file?(f ilename) * 
writoClfl, ’ (a32) ’ ) f oof ilo 
endif , 

if (runtypo) 100,200,300 

C EIGENHODE INPUTS 


100 print a, ’Busbar of modos desired: ’ 

if (runfile .aq. 'y') road 20, comaant 
road *, neig 

if (runfila .aq. *y') read 20, comment 
print a, 'Output File Hama:’ 
read 20, outfile 
open(13,file a outfile) 
if (savin .aq. ’y') than 

vrita(16, ’(a36)’) ’a Humber of modes desired? 1 
vrite(l6,’(i4)’) naig 

vrite(16,’(a33)’) ’a Output file?(f ilename) ’ 
writ a (16 , ’ (&32) ’ ) ontfila 
andif 


call READ FEB 


goto 1000 
C CSI INPUTS 

200 print a ,» Controller Definition File Name:' 
if (runfile .aq. ’y') read 20, comment 
read 20, confile 
open(l2 ,file s conf ile) 

print a, ’Please input type of control:’ 
print * 

print *, ’ 1. Full State Feedback' 

print a, » 2. Luenbarger Observer (1.1=0) ’ 

print a, » 3. Kalman Filter' 

print * 

if (runfila .aq. *y’) read 20, comment 
read *, contype 
contype a contype - 2 

print *, ’Initial time, final time, control-on time, step size: 
if (runfile .eq. ’y’) read 20, comment 
read a, tO,tf,tc,h 

print a, 'Forcing function ID, scale factor, damping coeff- a.b 

if (runfila .aq. *y*) read 20, comment 

read «, forceid.scalef ,adamp,bdamp 

print *, ’Phase lag fix?(y or n):' 

if (runfile .aq. ’y’) read 20, comment 
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read 21, tamp 

if (tastp ,eq. 'y' ) nolag * .true, 
if (temp . eq. 's') nolag ■ .false, 
print a, 'Gaia scale factors (4 total):* 
if (ruafile .eq. *y’) read 20, cement 
read * , qalph* 0 qbeta , qalphao , qbetao 
if (savin .eq. *y’) then 

write(16, *(a42)’) ** Costroller file name?(f ilename) * 
write(16,’(a32)*) eosfils 

write(16,'(a42)') '• Please iapat type of costrol: * 

write(16, *(iiO)') costype * 2 

write(16,*(a46)’) ’* Xaitial.final.eontrol-oa.step size?* 
write( 16 , » (4f 14 . 8) * ) t0,tf,tc,h 

vritq(16,’(a49)*) ** Forcing fsactioa, scale f, damping a.b? 
write(16,*(i4,flS.6,2fl2.8) *) forceid,scalef ,adamp,bdamp 
write(16, ' (a32) * ) •* Phase lag fix?(y or a)’ 
if (aolag) write(16,'(al)*) *y* 
if (.aot. aolag) write(16, *(al)*) 'a' 
write(16,*(a40)') *« Gaia seals factors (4 total)?* 
writ e( 16 , ’ (4f 14 . 8) * ) qalpha , qbeta , qalphao , qbetao 
eadif 

call READFEH 
goto 999 

C TRANSIENT RESPONSE INPUTS 


300 print *, 'Initial time, final time, step size:' 
if (runfile .eq. *y*) read 20, comeat 
read *, tO.tf.b 

print *, 'Forcing fsactioa ID, scale factor, damping coaff- a.b 
if (raafile .eq. *y') read 20, comeat 
read «, f orceid , scalef ,adamp,bdamp 
if (savia .eq. *y*) then 

write(16, *(a48)') ’* Initial, final, step size?* 
write( 16 , * (3f 14 . 8) ’ ) tO.tf.h 

write(16,*(a49)*) ’* Forcing function, scale f, damping a.b? 
write(16, *(i4,flS.6,2fl2.8) *) f orceid, scalef ,adamp,bdamp 
eadif 

call REAOFEM 


goto 999 

C OUTPUT OPTIONS 

999 print *, 'Output File Name:* 

if (ruafile .eq. *y*) read 20, comeat 
read 20, ostfile 
open( 13 , f ile»outf ile ) 

print *, ’Number of displacement results to output (max 10):' 
if (ruafile .eq. *y') read 20, comeat 
read *,ndisout 
do 600 i s l,ndisout 

print * , ’ Input node *, dot tor displacement output#' ,i 
read *,outlabel(i),outlabel(i+10) 

500 continue 

print *, 'Number of velocity results to output (mar 10):' 
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if (runfila .aq. * jr' ) road 20, comment 
read * ,nvelout 
do 600 iai.nvelout 

print * , ' Input nod* #, dof for velocity outputs', i 
road * ,outlabel(i+20) ,outlabel(i+30) 

600 eontinuo 

print *,'Send output ovory boo aany stops?' 
if (ruafilo .aq. *y') road 20, comment 
road * , outskip 
if (savia .oq. ’y‘) tboa 

vrito(16, '(a38)') '* Output fils aaao?(filoaaiBo) ' 
write(16,'(a32)*) outfilo 

writs (16, ’(a42) *) '» Humber of displacement outputs?' 
vrite(16, *(i4)*) ndisout 
do 660 i*l .adisout 

write(16,'(2i8)') outlabal(i),outlab*l(i*10) 

660 eontinuo 

srito(16, ’(a38)') '* Huabor of volocity outputs?' 
mrite(16, *0L4)O avolout 
do 660 ial,nvelout 

writo(16, * (2i8) ' ) outlab*l(i*20) , outlabal(i*30) 

660 continue 

orito(16, *(a44)*) ’* Send output ovary bow many stepj?’ 
srite(16, ’ (i3) ’) outskip 
ondif 

C AHIHATIOH 0PTI0H 

print * , * Animation Output? (y or n) : ’ 
if (runfila .aq. 'y') road 20, commant 
road 21, tamp 

if (tamp .aq. *y') animats « .true, 
if (tamp .oq. 'n') animats » .falsa, 
if (animats) tban 

print*, 'Animation fils name (filename) ' 
if (runfila .aq. 7 'y') road 20, comment 
road 20, animfilo 
open (15, file*animfilo) 
andif 

if (savin .aq. *y') tboa 

vrite(16,*(a41)*) ’• Sand animation output?(y or n)' 
if (animats) vrito(16,'(al)') 'y' 
if (.not. animats) vrito(l6,’(al)') 'n' 
if (animats) tboa 

vrito(16, ' (a31) ') ’* Animation file name?' 
vrito(16,'(a32)*) animfilo 
andif 
andif 

delta * b/2. 
dolsq a dolta**2 

1000 return 
end 


File: pmvmul.f 
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CsModale PMVHUL 
CsAntkor g. Airis, 
CsDate May 1990 
C»Bleek Fortran 


C 

C Subroutine PKVHUL 
C 

C Purpose: 

C T&ie subroutine multiplies a matrix ia vector lorn 

0 aad a vector. 

C 
C 
C 
C 
C 
C 
e 
c 
c 
c 


Arguments 

a - matrix ia vector form 

b - vactor 

c - reaslt vactor 

aeq - order of vactor aad square matrix 

jdiag - array of diagonal addresses for a 


C aabroatiaa PHVHDLCa, jdiag, b.neq.c) 


recursive aabroatiaa PHVH0L(a, jdiag, b, aaq, c) 

raal*8 a(l), b(l), c(l) 
integer jdiag(l), aaq 

do 100 i®l ,aaq 

c(i) » a(jdiag(i))*b(i) 

100 coatiaaa 


do 200 ia2,aaq 

do 300 jajdiag(i-l)+l,jdiag(i)-l 
k » jdiag(i) - j 
c(i) ■ c(i) ♦ a(j)*b(i-k) 

300 coatiaaa 

200 coatiaaa 

do 250 ia2,aaq 

do 400 jajdiag(i**l)+l , jdiag(i)-l 
k ■ jdiag(i) - j 
c(i-k) ■ c(i-k) ♦ a(j)ab(i) 

400 coatiaaa 

250 coatiaaa 

ratora 

aad 

C 

C 

C 

C Subroutiaa PMVHAD 

C 

C Porpoaa : 

C Multiply a. matrix ia vactor form aad a vactor and add the 

C. raaultaat vactor multiplied by a constant to a second vactor 
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multiplied by a second vector 


} 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


•100 


300 

200 


400 

250 


Arguments 

a 

b - 
e 

factl - 
fact 2 - 

jdiag - 
neq 


matrix ia vactor form 

vactor to ba multiplied with matrix 

rasaltaat aad vactor to ba added 

constant multiplier of matrix and first vactor 

constant multiplier of sacond vactor 

array of diagonal addresses for matrix 

order of vectors aad matrix 


subroutine PMVMAD(a, jdiag, b.neq, factl ,c,fact2) 

recursive subroutine PMVHAD(a, jdiag, b, neq, factl, c,fact2) 

xeal«8 «&)„ b£l},c&), fasti ,fact2 
integer jdiag(l), neq 

do 100 i«l,neq 

c(i) ■ fact2«c(i) ♦ factl*a(jdiag(i))*b(i) 
continue 

do 200 1*2, neq 

do 300 j»jdiag(i-l)*l, jdiag(i)-l 
k - jdiag(i) - j 
c(i) a c(i) ♦ factl*a(j)*b(i-k) 
continue 
continue 

do 250 ia2,neq 

do 400 j»jdiag(i-l)+l,jdiag<i)-l 
k a jdiag(i) - j 

c(i-k) » c(i-k) ♦ factl*a( j)*b(i) 
continue 
continue 

return 

end 


\ 


File: prepfem.f 


CaHodule PREPFEH 

C=Purpose Preprocess Structure Finite Element module for ACSIS 
C= Author K. Alvin 
CaDate May 1990 
C=Block Fortran 

C 

C 

C Subroutine PREPFEH 
C 

C Purpose: 

C This subroutine prepares the finite element mass, 

C stiffness, and S matrices in reduced profile vector form 
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c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 


Local variable*: 

sk Elaaaat Stiffness aatrir 

an Element Xui Matrix 

3a Local/Global DOF Mapping vector 

nseq Humber of olaaeat degrees of freedom 

ea,ep Material aad Propety id * for eleaeat 

subroutine PREPFE* 


include ’ shared. inc' 


C LOCAL VARIABLES 

paraaetor (MAXSEQ*24) 

ml *8 akOUlSEQ.MAZSEgj.saauZSSq.JUJSSQ^^e^e 
integer la(MAXSEQ) ,nseq,em,ep 

call REBUM 

C Set up akyliaa storage profile for global matrices 

call PR0FXLE(ix,id,jdiag,nnp,ael,4,6,mlen,ndof .mask) 

C Perform automatic domain decomposition 
call DOMDEC 

C Check six* of skyline profile against storage limitation 
if (alen .gt. MAXMLEB) then 

print*, 'PREPFEM: error, global matrix exceeded max. size' 
endif 

C Zero Global Matrices prior to assembly 

call ZEROVECT(stif ,mlen) 
call ZEROVECT(maas .mien) 

C ASSEMBLE EACH ELEMEBT MASS AMD STIFFBESS 

do 100 n*l,nel 

do 20 k*l,4 
j*ix(k,n) 

if ((etype(n).eq.l).and.(k.gt.2)) j » 0 
do 30 i*l,6 
kk=6*(k-l) ♦ i 
if (j .ne. 0) then 
la(kk) * id(i,j) 
else 

la(kk) * 0 
endif 

30 continue 

20 continue 
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if (atypa(a) . aq. 1) tkaa 
aaaq ■ 12 
•a ■ aaat(a) 

«p * aprop(a) 

call BEAM3D(n,ix(l,a) ,ix(2,a) t ix(3,n} ,coxyx,mat(l,8m) , 
aat(2, an) ,mat(3,aa) , prop(i, ap) ,prop(S,ep) ,prop(6,ep) , 

. prop (2 , ap) ,prop(3, ap) ,prop(4 , ap) , pia( 1 ,a) , sk , an) 

alaaif (ix(l,a) .a«. 0) tkaa 

print*, 'PREPFEH.-Eleaaat t jrpa not found, a 3 ' , a, ^types’ ,atype(a) 
aadif 

C ADD ELEMENT TO GLOBAL MASS AMD STIFFNESS 

call ADDSTF(ak,lm,atif , jdiag.aaaq) 
call ADDSTF(an,la,naaa, jdiag.aaaq) 

C SAVE THE ELEHEHT STIFFNESS FOE E-8Y-E COMPUTATIONS 

call SAYESE(ak*a,JU«q) 

100 coatiaaa 

C ADD LUMPED INERTIAS TO GLOBAL MASS? 

do 12S i=l,nnp s 

do 130 j«l,fl 

it (id(j,i) -«q- 0) goto 130 
k»jdiag(id(j,i)) 

maas(k) > aaaa(k) ♦ iaartia(j.i) 

130 coatiaaa 

12S coatiaaa 

C ASSEMBLE AND FACTORIZE aa (S MATRIX) 

me a 1. ♦ dalta«adamp 
kc » dalta*bdamp ♦ dalaq 
do 200 i 3 l,nlaa 

aa(i) 3 mc*maaa(i) ♦ kc*atif(i) 

200 coatiaaa 

♦ 

call SQLVER(aa,gn,jdiag,adof ,1) 

C INITIALIZE DISPLACEMENT AND VELOCITY VECTORS 

do 300 i 3 1 , aap 
do 3S0 j 3 1,6 

if (id(j,i) -aa. 0) thaa 
q(id( j ,i)) 3 qO(j,i) 
qdot(id( j ,i)) 3 qdotO(j.i) 
aadif 

3S0 coatiaaa 

300 coatiaaa 

ratora 

aad 

aabroatina SAVESX(ak,a,aaaq) 
iaclada ’ akarad.iac’ 
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real*8 sk(aseq,i) 
integer a.aseq 

k*0 

do 10 j»i„aseq 
do 20 1*1, j 
k»kel 

*eti£®OE,a)*skd,j) 

20 eoatiaae 

10 co&tiaae 

return 

•ad 

subroutine DOHDEC 

iaclado ’ shared. iac' 

logical aehk,adehk(HAXHODE,)UXDOK) 
integer ados 

do 10 j-l.lttXDOH 
a«ld(j)«0 
do 20 i*i,aap 

adebkd , j)» .false . 

20 eoatiaae 

10 eoatiaae 

ados&iaeO 

do 100 aal.nel 

ndo®*0 

nchJc»0 

do 200 vhile (achk.eq.O) 
ndom»ndom+l 

if (adoa.gt.adoaaia)- adoaaiaaadon 
aehk*l 

if (adchk(iz(l,a) .adorn/) nchk=0 
if (adcltic(ix(2,a),adom)) achk*0 
if (achk.eq.l) then 
•ldo»(n)*adom 
adcUcdzd.a) , adorn)*. true. 
adehk(ix(2,a) ,adon)».tnie. 
end if 

200 eoatiaae 

ueld(adon) *aeld (adorn) ♦! 
elaan(aeld(adoa) , adorn) *a 

100 eoatiaae 

retara 

ead 


File: profile.f 
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•» * 


C=Module PROFILE 

CsPurpoee Compute the number o f equations and set profile for K 
C»Autkor Bob Taylor 
C=Date wbo knove 

CsOpdate January 1939 by E. Pramono 
C= Block Fortran 


subroutine PROFILE ( ix,id, jdiag, nnp , nel , uen , adof , and , neq , mask) 

C+ ♦C 

C PURPOSE: C 

C THIS SUBROUT IHE COMPUTES THE HUMBER OF EQUATIONS REQUIRED C 

C TO SOLVE THE PROBLEM BT EL IH HUT INS RESTRAINED DEGREES OF C 

C FREEDOM FROM THE SYSTEM OF EQUATIOHS. KNOWING THE EQUATION C 

C NUMBERS COORESPOND IHG TO THE HODAL DEGREES OF FREEDOM, THE C 

C DIAGONAL ELEMENT LOCATIONS CAN BE COMPUTED FOR STORING THE C 

C GLOBAL STIFFNES MATRIX IH COMPACTED VECTOR FORM. C 


C 

C ARGUMENTS 
C 

integer ix(aea,l), id(ndof,l), jdiag(l) 
integer nnp, nel, nad, aeq, mask(l) 

C integer nnp, nel, nen, adof, nad, aeq, mask(l) 
C 

C .LOCAL ARGUMENTS 
C 

integer i, j. k. 1, a, n, jl, kl, 11, ml 
C 

C SET UP EQUATION NUMBERS 
C 

aeq a 0 

do 30 n > 1, nnp 
do 20 a a 1, adof 
j » id(a,mask(n)) 
if (j . eq. 1) goto 10 
neq » aeq + 1 
id(a,mask(n)) ■ neq 
jdiag(neq) a 0 
goto 20 

10 id(a,maak(n)) a o 

20 continue 

30 continue 
C 
C 

C COMPUTE COLUMN HEIGHTS 
C 

do 80 n a i, ael 
do 70 a ■ 1, nea 
ml a ix(m,a) 
if (ml .le. 0) goto 70 
do 60 1 ■ 1, adof 
11 a id(l,al) 
if (11 .eq. 0) goto 60 
do SO k ■ a, nea 
kl > ix(k,a) 
if (kl .le. 0) goto SO 
do 40 j si, adof 
jl » id(j.kl) 
if (jl .eq. 0) goto 40 
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i a HAXO(li,jl) 

jdiag(i) = MAXO( jdiagCi) , XABS(H“ji)) 

40 continue 

SO continue 

SO continue 

70 continue 

80 eeatiaue 

C 
C 

C COMPOTE DXAGOHAL PO UTTERS 
C 

aad * I 
jdiag(i) ■ 1 
if (neq . cq. 1) return 
do 90 i * 2, neq 

jdiag(n) » jdiag(a) ♦ jdiag(a-l) ♦ i 
90 continue 

aad a jdiag(neq) 

C 

return 

•ad 

C=Ead Fortran 


File: read.f 


CaRodule READ 
C= Author K. Alvin 
C^Date May 1990 
CaBlock Fortran 

C — 

c 

C Subroutine READFEM 

C 

C Purpoa* : 

C This aubroutin* roads fh« data file for the finite 

C element model. 

C 

C 

C Arguments 

C ctype - stores code for type of line 

C 

subroutine READFEM 
C GLOBALS 

include 'shared.iac' 

C LOCALS 

integer j ,a,ctype,GETTYPE 
chaxacter*132 aline 
real*8 in 

C INITIALIZE SIZE OF PROBLEM 
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sap a 0 
aal * 0 
ndof a 0 
adoaaia a 0 


C IDENTIFY CARD TYPE ABO ASSIGN DJPUT 


10 

100 

150 


raad(U,1000,aada9999) alias 
ctypa » GETTYPE( alias) 
if (ctypa) 10,10,150 


if (aliaa(l:4) .sq. 
if (aliaa(l:4) .sq. 
if (alias(l:4) .sq. 
if (alias(l:4) .sq. 
if (alias(l:4) .sq. 
if (alias(l:4) .sq. 
if (aliaa(l:4) .sq. 
if (alias(l:4) .sq. 
if (alia*(l-:4) .sq. 
if (alias(l:4) .sq. 


* HODS') goto 200 
’TOPO’) goto 300 
'ATTR') goto 400 
'MATE’) goto 500 
* PROP ’ ) goto 800 
’Fill’) goto 700 
'nrrr) goto soo 

•INER’) goto 900 . 
>EH9 >) goto 10 
’MESH’) goto 10 


priat a.'REAOFEH: Uaracogaizad cud typo; 
goto 10 


’ ,aline(l:4) 


C READ MODES 


200 road(ll,1000,sada9999) alias 

ctypa a GETTYPE(alias) 
if (ctypa) 200,250,100 

250 rsad(alias ,*) a, (coxyz(j ,a) , j»l ,3) 
if (a .gt. aap) nap a a 
goto 200 

C READ TOPOLOGY 

300 r aad( 11 , 1000 , aad»9999 ) alias 
ctypa » GETTYPE(aliaa) 
if (ctypa) 300,350,100 . 

350 rsad(aliaa,«) a.atypa(a) ,(ix(j ,a) , j»l ,4) 
if (a .gt. asl) aal a a 
goto 300 

C READ ATTRIBUTES 

400 raad(ll,1000,aada9999) alias 

ctypa ■ GETTYPE(aliaa) 
if (ctypa) 400,460,100 

450 raad(alias,«) a,aoat(a) , eprop(a) , (pia(j ,a) , j=l ,6) 
if (aldon(a).gt.adooaia) adoaaia s eldoa(a) 
goto 400 

C READ MATERIAL 

500 raad(ll,1000,aada9999) aliaa 

ctypa a GETTYPE(aliaa) 
if (ctypa) 600,550,100 

550 raad(aliaa,«) a,(n>at(j,a),jai,3) 
goto 500 

C READ PROPERTIES 
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600 r®ad(ll,1000,«nds9999) alia® 

ctyp® » GETTYPE( aline) 
if (etype) 600,660,100 
660 r®ad(aliac,*) a,(prop(j,a),j»l,6) 
g ota 600 

C BEAD Film 

700 r®ad(ll,1000,®ndo9999) alia* 

etype » OETTYPE(allae) 
if (etype) 700,760,100 
760 r®ad(alia«,*) n,(id(j ,u) , j«l ,6) 
goto 700 

C BEAD ramiL C0HDZTI0HS 

800 read(ll,1000,«nd<>9999) alia® 

etype • SEmPE<a2ia«) 
if (ctyp®) 800,860,100 
860 r«ad(alin®,*) n,j,q0(j,n),qdot0(j,a) 
goto 800 

C READ IBERTIA 

900 read(ll,1000,endo9999) alia® 

etype a GETTTPE(alia®) 
if (ctyp®) 900,960,100 
950 read(alia«,*) a.j.ia 

in«rtia( j ,n)»inertia( j ,a)*in 
goto 900 

1000 foraat(al32) 

9999 coatiau® 

return 

and 

C 

C 

C Snbroatia® READC0H 

C 

C Parpoa®: 

C This anbrontia® reads th® actuator and sensor 

C locations and th® gains for th® control system 

C 

C 

C Arguments 


c 

c 

ctyp® 

* 

stores cod® for type of line 

c 

c 

Abbreviations 

c 

HACT 

- 

nnmbor of actuators 

c 

HSEH 

- 

number of sensors 

c 

BKAT 

- 

locations of actuators 

c 

HD HA 

- 

array of displacement sensor locations 

c 

HVHA 

- 

array of velocity sensor locations 

c 

F1GA 

— 

control gain matrix 
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C F2GA ~ control gain matrix 

C L1GA - atata aatimator filter gain matrix 

C L2GA - atata aatimator Iiltar gain matrix 

anbrontina READCOH 

induda ’ ahared . inc ' 

C LOCALS 

raal*8 val 1 

intagar j ,n, ctypa, GETTYPE 
character* 132 alina , 

bval a 0 
hdval a 0 
hvval » 0 
hdbv&l » 0 
hvbval ■ 0 

C IDENTIFY CARO TYPE AND ASSIGN INPOT 

10 raad(12,1000,end*9999) alina 

100 ctypa » GETTYPE(alina) 

if (ctypa) 10,10,150 

150 if (aline(l : 4) .aq. 'NACT') goto 200 

if ( alina(l : 4) .aq. 1 NSEN 1 ) goto 300 
if (aline(l:4) .aq. 'BRAT') goto 400 
if (aline(l : 4) .aq. ’NONA*) goto 500 
if (alina(l:4) .aq. ’HVRA* ) goto 600 
if (alina(l:4) .aq. 'F1GA') goto 700 
if (alina(l:4) .aq. >F2GA>) goto 800 
if (alina(l:4) .aq. *L2GA*) goto 900 
if (alina(l:4) .aq. 'LiGA') goto 1100 
if (alina(l:4) .aq. 'END ') goto 10 

print «,'READC0N: Unrecognized card type; ' ,alina(l:4) 
goto 10 

C READ INPUT CARDS 

200 raad(12,1000,end;>9999) aline 

ctypa a GETTYPEC alina) 
if (ctypa) 200,250,100 
250 r«ad(aline,*) nact 
goto 200 

300 raad(12,1000,and a 9999) alina 

ctypa a GETTYPEC alina) 
if (ctypa) 300,350,100 
350 readCaline,*) naan 
goto 300 

400 raad(12,1000,and39999) alina 

ctypa a GETTYPE(alina) 
if (ctypa) 400,450,100 
460 read(aline,«) i,j,n,val 
bval 3 bval * 1 
b(bval) 3 vai 
brov(bval) a id(j,i) 
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bcol(bval) a a 
goto 400 

600 raad(12,1000,aad«9999) alia* 

ctypa ■ GETTTPE(aliaa) 

. if (ctypa) 500,660,100 
660 raad(alias,«) i,J,a,»al 
bdval a bdval ♦ 1 
bd(bdval) a v&l 
hd7o«(bdTal) a a 
bdeol(bdval) ■ ld(j,i) 
goto 600 

600 road(12 , 1000 , sad»9999) alia* 
ctypa ■ GETTYPEC alias) 
if (etypa) 600,680,100 
650 rsad(aliaa,«) i,j,a,val 
bwal ■ hwal ♦ 1 
hv(bvval) a -ral 
hvroo(hwal) » a 
bveol(brwal) a id(j,i) 
goto 600 

700 rsad(12,1000,sada9999) alia* 

ctypa a GETTYPEC alias) 
if (ctypa) 700,750,100 
760 raad(alias,*) i,j,a,val 

fl(a,id(j,i)) » qalpha*val 
goto 700 

800 rsad(12, 1000 . snd=9999) alias 
ctypa a GETTYPE(aliaa) 
if (ctypa) 800,860,100 
860 raad(aliaa,«) i,j,a,val 

f2(a,id(j ,i)) a qbata*val 
goto 800 

900 raad( 12 , 1000 , end a 9999 ) alias 
ctypa * GETTYPEC alias) 
if (ctypa) 900,960,100 
950 rsad(aliaa,«) i,j,a,val 

12(id(j ,i) ,a) a qbatao*val 
goto 900 

1100 raad(12,1000,sada9999) alias 

ctypa a GETTYPE(aliaa) 
if (ctypa) 1100,1150,100 
1150 raad(aliaa,«) i,j,n,val 

llCidCj ,i) ,a) a qalpbao«val 
goto 1100 

1000 fonaat(al32) 

9999 coatiaaa 
ratnra 
sad 

C 

C 
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C Function GETTYPE 

C 

C Purpose : 

C This function identifies whether a line is a character 

C input, data input or comment . 

C 

C 

C 

function GETTYPE(string) 

C GLOBALS 

character*132 string 
C LOCALS 

integer GETTYPE, ctype( 10) 
character*! head(lO) 

data head /•!*,*•», 't* ,'$* ,»X» , **»,»•* ,»C» ,'c* .» */ 
data etype ,-1,-1 ,-1,-1, -1,-1 ,0/ 

C LOGIC 

GETTYPE* 1 
do 100 i*l,10 

if (string(l:l) .eq. head(i)) GETTYPE=ctype(i) 

100 continue 

return 

end 


File: solver. f 


C*Kodule SOLVER 

C=Purpoae Solves the system of linear symmetric equations 

C= Author who knows 

C=Date 


C=0pdata January 1989 by E. Pramono 
C=Block Fortran 

C SUBROUTINE SOLVER ( BK , BR , JD I AG , HEQ , IFLAG ) 

recursive SUBROUTINE SOLVER ( BK , BR , JD I AG , NEQ , IFLAG ) 


C+ 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


PURPOSE: 

THIS SUBROUTINE SOLVES THE SYSTEM OF LINEAR SYMMETRIC 
EQUATIONS IN VECTOR FORM USING THE CROUT REDUCTION 
METHOD. 

ARGUMENTS: 

BK - GLOBAL STIFFNESS EQUATIONS IN VECTOR FORM 

BR - GLOBAL LOAD VECTOR 

JDIAG - LOCATION VECTOR FOR DIAGONALS IN [BK] 

NEQ - NUMBER OF EQUATIONS 

IFLAG - FLAG INDICATING WHICH FUNCTION IS TO BE PERFORMED 

1 -> FORWARD REDUCTION 

2 -> BACKWARD SUBSTITUTION 


+C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 


C 

c 
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Q ARGUMENTS 


BI<1), BE(1) 

INTEGER JDIAG(l), NEQ, I FLAG 
C LOCAL VARIABLES 

REAL*8 ZERO, EZERO, TOL, DAVAL, DOT, D, HDD, DO 

integer ldflag, jr, j, jd. jh, is, ie, k, jdt 

INTEGER JJ, ID, I, IB, IB 


C 

C 

C 

C 


c 

c 

c 

c 


c 

10 


c 

c 

c 

c 


c 

c 

c 

c 


20 

c 

c 

c 


JJ a 6 


NEW PARAMETERS 


2ERS • 0.000 
EZERO > 0 . 3D- 14 
TOL * O.SD-7 
LDFLAG » 0 


FACTOR BK TO UT*D*U OR REDUCE R 


JR a 0 

do 70 j ■ l, req 
JD a JDIAG(J) 

JH « JD - JR 

IS a J - JH ♦ 2 

IF (JH - 2) 60, 30, 10 

IF (IFLAG .HE. 1) GOTO GO 
IE a J - 1 
K a JR ♦ fi 
ID » JDIAG(IS-l) 


IF DIAGONAL IS ZERO COMPUTE A HORM FOR SINGULARITY TEST 


JDT ■ JDIAG(IE) ♦ 1 

IF (BK(JD) .EO. ZERO .AND. IFLAG .EQ. 1) THEN 
CALL DATEST (BK(JDT), JH-2, DAVAL) 

END IF 


REDUCE ALL EQUATIONS EXCEPT FIRST ROW AND DIAGONAL 


DO 20 I a IS, IE 
IR « ID 
ID a JDIAG(I) 

IH « MIN0(ID-IR-1,I-IS*1) 

IF (IH .GT. 0) BK(K) » BK(K) - DOT(BK(K-IH) , BK(ID-IH), IH) 

K ■ K * 1 
CONTINUE 


REDUCE FIRST ROW AND DIAGONAL 
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o o o o noon noon 


C 

30 


40 


IF (IFLAG .HE. 1) GOTO 50 
IE ■ JE ♦ 1 
IE » JO - 1 
K ■ J - JO 
DO a BK(JD) 

DO 40 I ■ IE, IE 
ID ■ JDIAG(K«T) 

IF (BKCID) .EQ. 0.0) GOTO 40 
0 » BK(I) 

BK(I) « BK(I)/BK(ID) 

BK(JD) ■ BK(JD) - D*BK(I) 
CONTINUE 


CHECK FOE POSSIBLE EBEOES AND PRINT WARNINGS 


ROD a BK(JD) 

IF (DABS (EDO) -LX. XOL*DABS(DB)) 

IF (DO .LT. ZERO .AND. ROD .GT. ZERO) 
IF (DO .GT. ZERO .AHD . ROD .LT. ZERO) 
IF (DABS(RDD) .LT. EZERQ) 


WRITE (JJ.2000) J 
WRITE (JJ ,2001) J 
WRITE (JJ ,2001) J 
WRITE (JJ ,2002) J 


COMPLETE RANK TEST FOR A ZERO DIAGONAL TEST 


IF (DD .EQ. ZERO .AHD. JH .GT. 0) THEM 

IF (DABS(RDD) .LT. TOL*DAVAL) WRITE (JJ.2003) J 

END IF 


REDUCE RIGHT HAND SIDE 

50 IF (IFLAG .Eq. 2) BR(J) » BR(J) - DOT(BK(JR*l) , BR(IS-l), JH-1) 

60 JR a JD 

70 CONTINUE 

IF (IFLAG .HE. 2) RETURN. 

C 

C DIVIDE BY DIAGONAL TERMS 

00 80 I a i ( HEq 
ID a JDIAG(I) 

IF (BK(ID) .HE. 0.0) BR(I) a BR(I)/BK(ID) 

IF (BR(I) .HE. ZERO) LDFLAG a 1 
80 CONTINUE 

C 

C 

C CHECK FOE ZERO LOAD VECTOR 

IF (LDFLAG .EQ. 0) WRITE(JJ,2004) 

C 

C BACK SUBSTITUTE 4 

J a HEQ 
JD a JDIAG(J) 

90 Da BR(J) 

J a J - 1 
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IF (J .LE. 0) RETURN 
JR > JDXAG(J) 

IF (JO - JR .LE. 1) GOTO 110 
XS»J-JD*JR*2 
I ■ JR - IS ♦ 1 
DO 100 la IS. J 

BR(I) a BR(I) - BK(I*K)*D 
100 COHTXHUE 
110 JO a JR 
GOTO 90 
C 

C WARNING FORMATS 

2000 FORMATC/M ! WARNING ! ! 1 - IS SOLVER, LOSS OF AT LEAST 7 DIGITS' 

♦ /18X, 'Of REDOC Hi G DIAGONAL OF EQUATXQH; ' ,4X,I5) 

2001 FORMAT(/M! WARMING !! 2 - XH SOLVER, SIGN OF DIAGONAL CHANGED* 

♦ /18X, 'WHEN REDUCING EQUATION; ' ,181,15) 

2002 FORMAT (/' ! ! WARNING 1 J 3 - IN SOLVER. REDUCED DIAGONAL IS EERO' 

♦ /18X, 'FOR EQUATION; *,2SX, IS) 

2003 FORMATC/M! WARNING M 4 - IN SOLVER, RANK FAILURE FOR A ZERO' 

♦ /18X, 'UNREDUCED DIAGONAL IN EQUATION; ' ,7X, IS) 

2004 FORMATC/M! WARNING ! ! S - IN SOLVER, ZERO LOAD VECTOR') 

C 

END 

C=End Fortran 
C=Modul« DATEST 
CaBlock Fortran 

C SUBROUTINE DATEST(A, JH.DAVAL) 

recursive SUBROUTINE DATEST(A, JH.DAVAL) 


c+ 



*C 

c 



C 

c 

TEST FOR RANK 


C 

c 



C 

c 

INPUTS; 


C 

c 

A(J) 

- COLUMN OF UNREDUCED ELEMENTS IN ARRAY 

C 

c 

JH 

- NUMBER. OF ELEMENTS IN COLUMN 

C 

c 



C 

c 

OUTPUTS; 


C 

c 

DAVAL 

- SUM OF ABSOLUTE VALUES 

C 

c 

c+ 



C 

+C 


C 

C ARGUMENTS 

REAL*8 A(l) , DAVAL 
INTEGER JH 

C LOCAL ARGUMENTS 

INTEGER J 
C 

DAVAL > O.ODO 
DO 10 J a i, JH « 

DAVAL a DAVAL ♦ DABS(A(J)) 
10 CONTINUE 
C 

RETURN 

END 
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c 

c 

Cs£od Fortran 

OKodulc 007 

C*81ock Fortran 
C FUNCTION D0T(A,B,H) 

recursive FUHCTIOH DOT(A.B.H) 


C* r — +C 

C PURPOSE: C 
C THIS FUHCTIOH SUBROUTXHE PERFORMS THE DOT PRODUCT OF TWO C 
C VECTORS. C 
C C 
C ARGUMENTS: C 
C A - FIRST VECTOR IHVOLVED IH DOT PRODUCT C 
C B SECOHD VECTOR IHVOLVED IH DOT PRODUCT C 
C H HUMBER OF ELEMENTS IH EACH OF THE TWO VECTORS C 


am*8 dot, r(a) 

INTEGER H 


DOT ■ 0.0 
DO 10 I » 1, H 

DOT ■ DOT ♦ A(I)*B(X) 
10 COHTIHUE 
C 

RETURN 

END 

CaEnd Fortran 


File: nophlag.f 


CaModalo HOPHLAG 
CaAuthor K. Alvin 
CaOata July 1990 
CsBlock Fortran 

C 

C Subrout in* HOPHLAG 

C 

C Purpoa* : 

C This subroutine solves for the structural displacement 

C and valocity vactora at tha hali-step lor the phase lag 

C correction loop, and gate nev maasuraaant 

C 

C 

C Arguments 

C dalbata - dalta * bdamp 

C 

subroutine HOPHLAG (zp) 

raal *8 zp(l) 
include ’ shared . inc ’ 

77 





c 


LOCAL VARIABLES 


iatagar i 

saal*« v(MAXDOF) .dalbota 
C LOGIC 

C ABB APPLIES FORCES TO RHS AMD PREPARE MASS MULTIPLIER 

da tO i>t,adof 

gs(i) ■ gs(i) ♦ f (i) 

v(i) » (t. * dalta*ada«p ) *q( i ) * dalta«qdot(i) 

10 eoatiaaa 

C SOLVE FOR RIGHT HARD SIDE, ga 

do 77 i » i,adoi 

gs(i) * d«lsq*ga(i) ♦ v(i)*masa( jdiag(i)) 

77 cm&xsatM 

If (bdamp .n« . 0.) t&aa 
dalbata ■ dalta«bdasp 

call PHVKAO(stif , jdiag.q.ndof , dalbata, gs , 1 .dO) 

•adif 

C SOLVE FOR DISPLACEMENT, q, OS IMG RHS AMD MATRIX S 

f 

call SOLVER(as,ga,jdiag,adof ,2) 

do 100 i*l,ndof 

»(i) » (ga(i)-q(i))/dalta 
100 coatiaaa 

call ZEROVECT(zp,aa«a) 

do 200 jj » l.hdval 
i » hdrow(jj) 
j * bdcol(jj) 

zp(i) « zp(i) ♦ bd(jj)*ga(j) 

200 coatiaaa 

do 2S0 jj « l.hmral 
i » baroa(jj) 
j » hacolCjj) 

zp(i) ■ zp(i) ♦ bv(jj)*v(i) 

2S0 coatiaaa 

ratara 

aad » 


File: . zerovect . f 


CsModala ZEROVECT 

CaParposa Iaitializa vactor of givaa loagtb to zaro 
Calathor X. Alvia 
CaQata May 1990 
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* r 


C 

C Subroutine ZEROVECT 

C 



C 

subroutine ZEROVECT(v.n) 

real*8 v(l) 
integer n 

do 100 isl ( n 
r(i) a o.dO 
100 continue 

return 

end 


File: lu.f 


SUB ROOT lire LUFACT(A,N, PIVOT .DET.IER, UMAX) 
0****«*e**«****«****«»«**M****«****»****»**«*«»«*******«»*«*»*< 

C 

C SOBROOTIHE FACTOR USES GAOS SI AH ELIMIHATION WITH 

C PARTIAL PIVOTIHG AHD IMPLICIT SCALING TO DETERMINE 

C THE L*U DEC0MP0SI0H OF A SQUARE HATRIX "A" OF 

C ORDER N. THE ALGORITHM ALSO FINDS THE DETERMINENT 

C OF "A". UPON COMPLETION, THE ELEMENTS OF THE UPPER 

C TRIANGULAR MATRIX "U" ARE CONTAINED IN THEIR RESPECTIVE 

C LOCATIONS IN MATRIX "A". THE ELEMENTS OF HATRIX "L" 

C ARE CONTAINED IN THE LOVER TRIANGULAR PORTION OF "A", 

C BUT ARE SCRAMBLED WITH RESPECT TO "U" BECAUSE OF ROW 

C INTERCHANGE OPERATIONS NOT PERFORMED ON THE ELEMENTS 

C OF "L". THE VECTOR PIVOT (SEE BELOV) MUST BE USED TO 

C UNSCRAMBLE "L" IF IT IS TO BE USED FOR OTHER OPERATIONS. 

C 

C VARIABLES: 

C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 

C*««*****«**«**e*e 

INTEGER PIVOT(l) ,IER,N,I, J,K,IO,NMAX 
REAL*8 A(NHAX,1) ,S(1000) ,C(1000) .DET.TEMP 
DET«1.0D0 

CM***W«SM*M**M»**M»*»**m»**«»***M**«**»*****MWM***M* 

C 

C FIND THE ROV NORMALIZING COEFFICIENTS S(I) FOR IMPLICIT SCALING. 

C EXIT ROUTINE IF ANY S(I)=0.0 


A=FULL SQUARE MATRIX (DOUBLE PRECISION) 

REORDER OF MATRIX A (INTEGER) 

PIVOT= VECTOR CONTAINING A RECORD OF 
ROV INTERCHANGES. THE INTEGER 
VALUE PIVOT(X) IS THE ROV VHICH 
VAS INTERCHANGED VITH ROV X AT 
FORVARD ELIMINATION STEP K. (INTEGER) 
DET»DETERHINENT OF MATRIX A (DOUBLE PRECISION) 
IER=ERROR FUG. IF IER=1, THE MATRIX A VAS FOUND 

TO BE SINGULAR, AND THE ROUTINE VAS EXITED. . IF 
IER=0, THE DECOMPOSITION VAS SUCCESSFUL. 
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nano 


C*®*®®®®*****^,* ®®«®®*®*®«*»tt®®*«®e®«*«®c»**«®a®s®®***®*®*®<!ii!<sa»*:S! 

DO 100 I»i,H 
S(I)»O.ODO 
DO 110 J*1,N 

IF (ABS(A(I,J)) .GT.S(I)) S(I)«ABS<A(I,J)) 

110 CONTINUE 

IF (S(I).Eq.O) THEH 
m«i 
DITsO.ODO 

asms 

END IF 
100 CONTINUE 

C®®**®®**®*®*®®*®*®*®®®®®®®*®®**®®**®®®*®**®*®**®®**®**®®*****®®* 


START FORWARD ELIMINATION STEP X 


i*®***«****®®®«****®® 

DO 130 X«l,Jhl 

C*®®®**®®®®*®®®®®*®®®*®®®**®*®®***®®®*®®*®*®®****®®*®®®®*®®®®®®®® 

c 

C DETERMINE PIVOT ELEMENT A(IO,K) BY FINDING THE ROW 10 

C BETWEEN K USD N CONTAINING THE MAXIMUM NORMALIZED . 

C VALUE IN COLUMN X. SET PIVOT(X)®IO 

C 

C®*®®®®®*®******®®®®®**®®®®®®***®®***®*®**®*®®**®*®®®®®®®*®®®®®®® 

C(K)«O.ODO 
DO 130 I*K,H 

TEMP«ABS(A(I,K)/S(I)) 

IF (TEMP. GT. C(lO) THEH 
C(K)»TEMP 
IO»I 
END IF 
130 CONTINUE 
PIVOT(K)®IO 

C®®*®®®®®®®®®®®®*®®®®®®®®*®®®®®®®®®*®**®®®*®®*®®®®®®®®®®®®*®®*®*® 

C 

C EXIT ROUTINE IF ALL VALUES IN COLUMN K AT OR BELOW 

C THE MAIN DIAGONAL ARE EqUAL TO 0.0 

C 

C*®®**®***®**®®®®®*®®*®®*®*®*****®®*®®®®®*®**®®®**®®*®*®®®®****®* 

IF (C(X).EQ.O.O) THEH 
IER»1 
DETaO.ODO 
RETURN 
END IF 

C**®**®®®®®**®®*®*®®®®®®*®®®®*****®®®®®®*®®*®*®*®®®®®®®®*®®®®®®*® 

c 

C INTERCHANGE ROWS 10 AND X FOR COLUMNS X TO N. SKIP IF IO=K. 

C SET DET«-DET IF ROWS ARE INTERCHANGED. 

C 

C**®**®*®***®®®*®*®®*®®®*®®®**®®®®®***®®*®®*®®®®***®®®*®*®®®®*®*® 

IF (lO.Eq.K) GOTO 150 
DET®-1.0D0*DET 
DO 140 J»K,N 
TEMP®A(K, J) 

A(K, J)®A(IO, J) 

A(IO,J)=TEHP 
140 CONTINUE 
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C*******************«*«*««******************v***«***!«e4*i3** «««:«** 

c 

C ELIMINATE counts X BELOV MAIN DIAGONAL BY MULTIPLYING 
C ROW K PROM COLUMH K TO H BY A(I ,K)/A(K ,K) AND SUBTRACTING 

C FROM ROW I. STORE THE MULTIPLIER FOR ROW I IN THE ELIMINATED 

C COLUMI X. MULTIPLY THE RUNNING PRODUCT DET BY DIAGONAL ELEMENT A(K,K). 

C 

C*************«*******«*a**»*****************************'**** , * , «** 


ISO DO 160 lag+l.N 

A(I,K)«A(I,K)/A(K,K) 

DO 170 J»K*1,N 

A(I,J)oA(I,J)-A(I,K)*A(K,J) 

170 CONTINUE 

160 CONTINUE 

DET»DET*A(K,K) 

120 CONTINUE 

C***********************.**************************************** 

c 

C CHECK LAST ROB/COLUMN FOR SINGULARITY. IF THERE IS NO ERROR, 

C COMPLETE CALCULATION OF THE DETERMINES! , SET THE ERROR FLAG 

C TO INDICATE NORMAL COMPLETION, AND EXIT. 

C 

C*****«***»*«******«*«*********«******************************«** 

IF (A(N.N).Eq.O.O) THEN 
IER»1 
DET=O.ODO 
RETURN 
END IF 

DET»DET*A(H,N) 

IER=0 

RETURN 

END 

SUBROUTINE LUSOLVE(A,N,B, PIVOT, NMAX) 
C**************************************************************** 

c 

C SUBROUTINE SOLVE 

C 

C**************************************************************** 
INTEGER PIVOT(l) ,N,I, J,K,NHAX 
REAL*8 A(NMAX,1), B(l), TEMP 
DO 100 Kal.N-t 

IF (PIVOT(K).EQ.K) GOTO 110 
TEHPaB(K) 

B(K)»B(J) 

B(J)aTEMP 

110 DO 120 I=K*1,N 

BCI)-B(I)-A(I,K)*B(K) . 

120 CONTINUE 

100 CONTINUE 

B(N)aB(N)/A(H,N) 

DO 130 IoN-1,1,-1 
DO 140 J»I*1,N 

B(I)»B(I)-A(I,J)*B(J) 

140 CONTINUE 

B(I)=B(I)/A(I,I) 

130 CONTINUE 
RETURN 
END 
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CsEHD Foams 
csBECI FACTA 

esPOaPQSg - Factors the A aatriz u L C « A, vita partial pivoting 
e*ASTSBE 8 X BELT IX, Sopt. 24 0 1987 
e 

e Input 

c aaat--— [ja X n] aatriz to bo factored, dastroyad on output 

e up — .-x—probla* siso 

e 

e Output 

e gnat—- contains tha LU dacoaposition 

c 

subroutino FACTA (anat ,np,arov,lp) 

c 

real*8 usat(*),eta 
tat agar lp(*) 

c 

do SO i«l,np 
SO lp(i)«i 
c 

c Find largest pivot . 
c 

do 100 k»l,np-l 
aaazkaO. 

BBSaxaO 

do 200 mek.np 

if (abs(aBat(lp(a)+(k-l)*nrou)) .gt. aaaxk) then 
aaaxks&bs ( anat ( lp(a) +(k-l ) *nrow) ) 


maz*s 

aadif 

200 continue 
l-lp(k) 

lp(k)alp(maaz) 

lp(anax)*l 

do 400 i=k+l ,np 

ata»anat ( lp( i ) ♦ (k- 1 ) *nrov ) /amat ( lp (k) ♦ (k- 1 ) »nr ow ) 
aaat ( lp( i ) •♦■ (k- 1 ) enrou) ■ et a 
do S00 jak+l.np 

anat(lp(i)*(j-l)*arou)saBat(lp(i)+(j-l)*nrou)- 
eta*amat(lp(k)+( j~l)*nrov) 

500 continue 

400 continue 

100 continue 


return 

and 


C EHD FORTRAB 
c=DECK LUSOLV 

c=PTJRPOSE - Solve L U x • b, 
c=AUTH0R 8 X' BELYIS, Sept. 24, 1987 
c * 
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c First solvss L j » b, than 0 x « y 
c 

c Input 

e mt*— [i X a] matrix factored by FACTA 

c np - — — problem six* 

c lp —pointer vector based on pivoting 

e rhs— RHS of aquation 

c 

c Output 

c aaat contains tha LU decomposition . 

c x— — — " the solution vector 

c 

subroutine LUSOLV(aaat ,np ,nrow , lp , rhs ,x) 
c 

raal*8 anat(e) ,x(*) ,rhs(») . 
integer lp(«) 
c 

do 54 i»l ,np 
x(i)<>rhs(i) 

50 continue 

c A - 

c Solve Lover System— — — — — — — — — 

c**a* Outer loop 

do 100 k*l,np 
rhs 00 =x (lp 00) 
if (rhs(k) . eq. 0.) go to 100 
c 

c**m Inner loop 
c 

do 200 i=k+l,np 

x(lp(i))=x(lp(i))-amat(lp(i)a(k-l)*nrov)*rhs(k) 
200 continue 
100 continue 
c 

c Solve Upper System— —————— — — — 

c 

do 300 k®np,l,-l 

x(k)=rhs(k)/amat(lp(k)+(k-l)anros) 
do 400 ial.k-1 

rhs(i)arhs(i)-x(k)*a«Bat(lp(i)+(k-l)*aroB) 

400 continue 

300 continue 
return 
end . 


File: prepcon.f 


C=Hodule PREPC0H 

C=PurpoBe Preprocess control analysis module for ACSIS 
C= Author K. Alvin 
C=Date June 1990 
CaBlock Fortran 

C 

C Subroutine PREPCOH 
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c 

C Purpose: ' 

6 This subroutine prepares ec, the control prediction integration 

C Matrix and eo, the observer construct matrix S (M+delta«‘D+delsq*K) 

C 

C 

C 

subroutine PREPCQS 

include ’ shared, inc' 

C LOCAL VARIABLES 

real* 8 me ,kc 
integer ier 

e Losie 

call READCQH 

C ! Form Control Prediction Integration Coefficient Matrix 
C 

C contype * -1 : Full State Feedback 

C contype » 0 Lnenberger Observer 

C contype a *1 Kalman Filter 

do 10 i ■ l.nact * naen 
do 20 j ■ l.nact ♦ nsen 
ec(i, j) ■ O.dO 
20 continue 

10 continue 

if (contype) 100,200,400 

J 00 ncsi * nact 

do 110 i ■ l.nact 
do 120 jj ■ l.bval 
j • bcol(jj) 
k » brov(jj) 

ec(i,j) » ec(i,j) ♦ delta*f2(i,k)*b(j j)/mass(jdiag(k)) 

120 continue 

110 continue 

goto 600 

200 ncsi * nact ♦ nsen 

do 210 i > l.nact 
do 220 jj » l.bval 
j * bcol(jj) 
k » brov(jj) 

•c(i, j) » ec(i.j) ♦ delta*f2(i ,k)*b( j j)/mass( jdiag(k)) 

220 continue 

do 240 j ■ aact+l.ncsi 
do 250 k > l.ndof 

ec(i.j) » ee(i.j) ♦ delta*f2(i,k)*12(k,j-nact) 
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250 

240 

210 


270 


290 

260 


400 


480 


S10 

800 

470 


830 


880 

520 

600 


1100 

C 


continue 

continue 

continue 

do 260 ii ■ l.hwal 
i ■ nact ♦ hvxow(ii) 
do 270 jj ■ l.bval 
j • bcolCjj) 
k » brou(jj) 

il (hvcol(ii) .ne. k) goto 270 

ec(i,j) ■ ec(i, j) ♦ delta*hv(ii)*b(jj)/naaa(jdiag(k)) 
continue 

do 290 j » nact*l,ncai 
k ■ kvcol(ii) 

ec(i,j) « ec(i, j) ♦ delta*h»<ii)*12(k, j-nact) 
continue 
continue 

goto 600 

ncsi s nact ♦ naen 

do 470 i 3 l.aact 
do 480 jj a l.bval 
j « bcol(jj) 
k » broe(jj) 

ec(i, j) a ec(i,j) ♦ deltaaf 2(i ,k)*b(j j)/maaa( jdiag(k)) 
continue 

do 800 j a nact+l,ncai 
do 510 k a i.ndof 

ec(i,j) a ec(i, j) ♦ delta*f2(i,k)»12(k, j-aact)/ 
mass( jdiag(k) ) 

continue 

continue 

continue 

do 520 ii a l.bwal 
i a nact + hvrou(ii) 
do 630 jj a i.bval 
j » bcol(jj) 
k » brow(jj) 

il (ivcol(ii) .ne. k) goto 530 

ec(i,j) » «c(i.j) ♦ dolta»hv(ii)*b( j j)/nass( jdiag(k)) 
continue 

do 850 j a nact>l,ncai 
k a bvcol(ii) 

ec(i,j) a ec(i, j) ♦ delta*hv(ii)*12(k, j-oact)/ 
naaa( jdiag(k)) 

continue 

continue 

continue 

do 1100 i * l.ncai 

ec(i,i) a «c(i,i) + l.dO 
continue 

FACTORIZE ec 

call FACTA(ec,ncai,KAXCSX .pivot) 
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if Cits . «q. 1) than 

prist « , ’ PREPCON : Singular Matrix for Control Integration* 
•ffidii 

C Fora 9b# error Integration Coefficient Matrix, eo 

k « 1. ♦ delta«adasp 
ke » delta«bdamp ♦ delsq 
do 1200 i»l.mlen 

ee(i) ® ac*aass(i) * kc*stif(i) 

1200 continne 

C FACTORIZE eo 

eall S0LYER(eo,go, jdiag.ndof ,1) 

C Initialize Observer States 

sail S8OTZCT (qe , sdef) 
eall ZEROVECT(qedot.ndof) 
eall ZEROVECT(pe.adof) 

return 

end 


File: control. f 


CsHodule CONTROL 
CsAnthor K. Alvin 
C=Date Hay 1990 
CaBlock Fortran 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


Subroutine CONTROL 
Purpose: 

This subroutine carries out the numeric integration of one time 
step of the control system. 


Arguments . 

q*p 

qedotp - 

PP 

z 


estimated displacement vector at half time step 
estimated velocity vector at half time step 
generalized momentum (f-D*qedotp~K*qep) 
measured sensor output 


subroutine CONTRQL(z) 


include * shared. inc* 
real*8 z(l) 

C LOCAL VARIABLES 


real*8 qep(HAXDOF) , qedotp (HAIDOF) .pp(HAXDOF) . v(HAIDOF) 
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C LOGIC 

C Fora BBS of Control Prediction Equation Sat 

C 

C contype ■ -1 Full State Feedback 

C contype ■ 0 : Lnenberger Observer uith Li = 0 

C contype ■ ♦! : Kalman Filter v/generalized momentum variable 

4 

call ZZBQVECKgc.ncsi) 

if (contype) 100,200,300 

100 continue 

do 110 i ■ l,ndof 

q*(i) ■ q(i) 

qedot(i) » qdot(i) 

110 continue 

200 , continue 

do 210 i ■ l.ndof 

qep(i) ■ qe(i) ♦ delta*qedot(i) 
qedotp(i) * qedot(i) 

pp(i) « f(i) - mass( jdiag(i))*adamp*qedotp(i) 
v(i) * qep(i) bdaop*qedotp(i) 

210 continue 

call PNVHA0(stif , jdiag, v.ndof ,-l .d0 ,pp,l .d0) 

do 220 isl.nact 
do 230 j » l.ndof 

gc(i) ■ gc(i) - fl(i,j)*qep(j) - f2(i, j)*(qedot( j) ♦ 
delta*pp ( j ) /mas e ( j diag ( j ) ) ) 

230 continue 

220 continue 

if (ncai . eq. nact) goto 600 

do 240 ianact+l.ncai 
k a i - nact 
gc(i) » z(k) 

240 continue 

do 245 ii » l,hdval 
i a bdrou(ii) ♦ nact 
j a hdcol(ii) 

gc(i) ■ gc(i) - bd(ii)*qep( j) 

24S continue 

do 250 ii a l'hwal 
i a kvrou(ii) ♦ nact 
j a hvcol(ii) 

gc(i) ■ gc(i) - hv(ii)*(qedot( j)>delta*pp( j)/masa( jdiag( j))) 
250 continue 

goto 600 

300 continue 

do 310 i a l.ndof 
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Vp(i) • q«(i) 

pp(i) » f (i) - aaBB(jdiagd))*adanp*qapd)/dalta 
▼(i) ■ (1 ♦ bdamp/dalta)*qapd) 

310 coatiaaa 

call P*7HJU>(«*if,jdiag,v,adof,-l.dO,pp,i.dO) 

do 320 i»l,aact 
do 330 j » l,adof 

gc(i) « gc(i) - 21(i,j)*q«p(j) - 22(i, j)*(pa(j) ♦ 

. dolta*pp(j))/Bass(jdiag(j}) 

330 coatiaaa 

320 eoatiaaa 

do 340 isaaet*l,acsi 
1 « i “ aaet 
ge(i) « *(k) 

340 coat lane 

do 346 ii ■ l.hdval 
i » bdxwdi) aact 
j « idcol(ii) 

gc(i) ■ ge(i) - bddi)*q«p(j) 

345 coatiauo 

do 350 ii ■ 1 ,hvval 
i « hvrop(ii) *■ aact 
j ■ hveol(ii) 

gc(i) ■ gc(i) - hv(ii)*(pa(j) ♦ dalta*pp(j))/ma.Ba(jdiag(j)) 
350 coatiaaa 

C FIHD r, C0HTR0L AHS STATE CORRECTION FORCES 

600 caU LOSOLV(ac,acai.HAXCSX. pivot, gc.r) 

do 610 j«l,aact 
a(j) • r< j) 

610 coatiaaa 

do 620 j*nact+l ,ncsi 
gaasa( j-aact) ■ r(j) 

620 coatiau* 

C FIHD CONTROL CONTRIBUTION TO RHS VECTOR FOR 
C OBSERVER AHD STROCTRE 

do 710 ial.adof 
gt(i) « O.dO 

god) ■ O.dO 

gk(i) * O.dO 

710 coatiaaa 

do 720 jjal.bval 
i » brow(jj) 
j » bcol(jj) 

gt(i) a gs(i) ♦ b(jj)*a<j) 

go(i) a go(i) ♦ b(jj)*a(j) 

gk(i) a gk(i) ♦ b(j j)*a(j) 

720 coatiaaa 

ii (coatjpa .aq. 0) than 
do 725 i > i.adof 

do 730 jai,asaa 

god) » god) ♦ maBB(jdiagd))*I2d,j)*gajnma.(j) 

730 coatiaaa 
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725 continue 

elseif C coat 7 p • .eq. 1) than 
do 73S i » l.ndof 
do 740 j»l,naen 

go(i) ■ go(i) + (12<i,j)anaaa(jdiagCi))aii(i,j) /delta) 
*gaaaa(j) 

gk<i) ■ gk(i) ♦ 12(i, j)*gaama(j) 

740 continue 

735 continue 

end if 

return 

end 


File: secorder.f 


CaModnle SECONDER 
C* Author K . Alvin 
ODate Nay 1990 
CsBlock Fortran 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Subroutine SECORDER 
Purpoae: 

Solvea the second-order dynamical equation: 

Nr" ♦ Ox' ♦ b • 1 ♦ j 

at time (a+1) given f(n*l/2), g(n+l/2) and x,x' at n by 
the midpoint implicit integration rule. 

Step size is 2*delta. 

0 is of the form (alpha«H ♦ beta*K), f is an applied force, 
and g is assumed to be other applied force from a feedback 
control loop. The matrix E is the factored fora of the 
integration coefficient matrix: E a (X ♦ delta*D ♦ delta'2*K) . 


Arguments : 

m 

k 

alpha 

beta 

f 

S 

e 

x 

xd 

delta 

delaq 

jdiag 

ndof 


matrix N 

matrix K 

scalar alpha 

scalar beta 

Force vector fCn+1/2) 

Feedback force vector g(n*l/2) 
martix S 

Variable vector x(n) 

Variable vector x' (n) 

Half of integration time step 
delta‘2 

Diagonal location pointer for X,K,E matrices 
Humber of equations and length of q 
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C t mass multiplier for RHS preparation 

C delbeta - delta * beta 

C subroutine SECORDER(a, k, alpha, beta, i,g,e,x,xd, 

C delta, delsq,jdiag,adof .KAX9QF) 

recursive subroutine SECORDER(a,k, alpha, beta.f ,g,e,x,xd, 
delta, delsq, jdiag.ndof .MAXDOF) 

C ARGUHESTS 

real *8 a(l),k(l), alpha, beta, f (1} ,g(l),e(l) 
real«8 x(l) ,xd(l) , delta, delsq 
integer jdiag(l) ,ndof .HAIDOF 

C LOCAL VARIABLES 

integer i 

real *8 ▼ (3000) ^delbeta 

C LOGIC 

C ADD APPLIED FORCES TO KHS ADD PREPARE MASS MULTIPLIER 

do 10 i»»l,ndoi 

g(i) ■ g(i) + f(i) 

v(i) ■ (1. ♦ delta«alpha)«x(i) ♦ delta*xd(i) 

10 continue 

C SOLVE FOR RIGHT HAHD SIDE, g 

do 77 i»l ,ndof 

g(i) a delsq*g(i) ♦ v(i)*m(jdiag(i)) 

77 continue 

if (beta .ne. 0.) then 
delbeta » delta*beta 

C Activate EBE computations lor internal force by using STIFFRC 
C subroutine. Otherwise use PKVHAD (profile matrix/vector mult-add 

C call PHVHAD(k, jdiag,x,ndof .delbeta, g, 1 .dO) 

call STIFFRC(x, delbeta, g) 
endif 

C SOLVE FOR DISPLACEHEHT, q, US IMG RHS AMD MATRIX E 
call SQLVER(e,g, jdiag.ndof ,2) 
do 100 i>l,ndof 

xd(i) ■ 2.*(g(i) - x(i))/delta - xd(i) 
x(i) ■ 2.*g(i) - x(i) 

100 continue 

return 0 

end 


File: measure. f 
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CaHodule MEASURE 
CaAuthor K. Alvin 
C=Data Kay 1990 
Csfllock Fortran 

C 

C 

C Subroutine MEASURE 

C 

C Purpose : 

C This subroutine storoa now measured aanaor data by using the 

C previous displacement and velocity vectors at the sensor 

C locations 

C 

C 

C Arguments 

C zp - measured sensor data array 

C 

subroutine KEASURE(zp) 

include ’ shared . inc ' 
real*8 zp(l) 

call ZERQVECT(zp,nsen) 

do 100 jj » l.hdval 
i » hdrow(jj) 
j ■ hdcolCj j) 

zp(i) ■ zp(i) ♦ hd(jj)*q(j) 

100 continue 

do 200 jj ■ l.hwal 
i » hvrow(jj) 
j * hvcol(jj) 

zp(i) » zpCi) ♦ hv(jj)*qdot(j) 

200 continue 

return 

end 


File: eigens.f 


C=Module EXCESS 

C=Purpose Find Eigenmodes given Mass, Stiffness Matrices 
C=Author K. Alvin 
CsDate March 1990 
C=Block Fortran 

C 

subroutine EXCESS 
C 

C 

C 

C COKHOS ABD GLOBALS 
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« 4 


iaeloda ’sharad.iac' 

C LOCAL VARIABLES 

para*atar(NVI*100 ,HHNVH*NAXDOF*HVM ,HXCHVsNVH*(NVM*l)/2) 
iatagar ialCNAXDOF^awae.i.j.k.kk^oat 
?«al*8 vl(NNNVN) ,vx(HNNVH) ,akk(NXCNV) t &ffl»(ilXCNY) 
saal*8 xx(NVH,NVl) ,eigv(NVH) ,aigold(NVH) 
real*8 tolaig,toljae, omega, fhz 

toleiga.OQOl 
owt * 13 

toljae ■ toleig 
avoe * mia0(2*aaig,100) 
awac « miaO(avec,adof ) 

C SET Iff ISL VECTOR 

iel(l) a i 
do 50 j«2,sdo£ 

i»Kj) a j - jdiag(j) ♦ jdiag(j-l) *1 
50 continue 

C CALL EIGENSOLVER 

call SSPACECstii.maaa.vl.irr.akk.aaffls.zr.eigv.aigold.isl, 
j diag , aaig ,av«c , ndof . tolaig , tol j ac , out ) 

C WRITE OUTPUT 

writaCoat.*) *EIGEH ANALYSIS RESULTS:* 

writaCout,*) • RADIAL CYCLIC* 

writ a (out, *) * MODE EIGENVALUE FREQUENCY FREQUENCY* 

writaCout ,*) * — — - — — --* 

writa(out,*) 
do 100 ial.aaig 

omaga a dsqrt(eigv(i) ) 

ihz a omaga/(2*3. 141692654) 

writaCout,’ (i5,3(3x,gl2.5))’) i,eigv(i) .omega, Ihz 
100 coatiaaa 

writ«(out,*) * EIGENVECTORS : ’ 
do 200 j*l,naig,5 
writaCout,*) 
do 300 i*l,ndoY 

k. « ndof*(j-l) ♦ i 

writaCout, *(i5,5(lx,gl2. 5))') i,(irr(kk) ,kk=k,k+4*ndof ,ndof) 
300 coatiaaa 

200 coatiaaa 

writaCoat,*) 'NASS MATRIX DIAGONAL:* ' 
do 400 i*l,nnp 
do 450 j*l,6 
if (id(j ,i) .ne.O) than 

writaCoat,*) i, j ,id(j ,i) ,mae«(jdiag(id(j ,i))) 
endif 

450 coatiaaa 
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400 


continue 


return 

end 

C=End Fortran 


File: singaig.f 


SUBROUTIHE SSPACE (AK,AM,VL,VR,AKK,AMM,XX,EIGV,EIGCLD ,ISL, 
1 IDIAG .HEIG.NVEC , HD0F , TOLEIG , TOLJ AC , HV) 


c 

e 

e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 

c 

c 

c 


input : 


AK 

AM 

TSL 


IDIAG 

MEIG ; 

HVEC 

NBOF 

TOLEIG 

TOLJAC 

HV 


stiffness matrix (profile values) (NDOF) 

consistent mass matrix (profile values) (HDOF) 
stores in position “i" the -row # of tip of 
column "i" (NBOF) 

position of diagonal terms in profile (HDOF) 

S of required eigenvalues 

9 of subspace vectors 

9 of degrees of freedom 

tolerance for eigenvalues convergence 

tolerance for Jacobi convergence 

logical unit number for output 


output : 

VL(HDOF.HVEC) 

- VL(.,1..MRH0D) : 

- VL(. .HRMQD. .HVEC) : 

VR(HDOF,HVEC) : 

- VR(.,1.. HRMQD) : 

- VR(.,HRM0D'. .HVEC) : 

AXX( HVEC* (HVEC ♦ l)/2) 
AMM(HVEC*(NVEC ♦ l)/2) 
XX(HVEC*HVEC) 

EIGV(HVEC) 

- EIGV(1..HRM0D) : 

- EIGV(HRM0D . .HVEC) : 


working array 
AM times rigid modes 
subspace at the previous step 

eigen-vectors 
rigid modes 

subspacs at this step (eigenvectors) 

stiffness matrix in the subspace 
consistent mass matrix in the subspace 
subspace eigenvectors 

current eigenvalues 
0 eigen-values 
following eigenvalues > 0 


EIGOLD(HVEC) : same as EZGV 


IMPLICIT REAL *8 (A-H.O-Z) 

DIHEHSIOH AK(1) ,AM(1) .VL(HDOF.HVEC) ,VR(NDOF,HVEC) ,AKK(1) , 
1 AMX(i) , II (HVEC* HVEC) ,EIGV(1) ,EIG0LD(1) ,ISL(1) , IDIAG(l) 

C 

VRITE (HV.1003) HEIG .HVEC, HDOF, TOLEIG 
C 

CALL IHVECT ( AK, AM, VL.VR, IDIAG, HDOF, HVEC) 

CALL FACT (AK, IDIAG, ISL.HDOF.HV) 

CALL HULL (AK, AM, VL.VR, IDIAG, ISL.HDOF.HRMOD) 
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c 

c HHHQD : S of rigid sodas 

C 

KRXTE (STB, 1004) HRMOD 
ffSUl*«VEC-IRIOD 

C - 

c m ORTHO (VR.YL.HDOF.HRSOD.HVEC) 

C 

SXtsO 

SSSAXsjg 

®VECl»lVEC-i 
DO S I»1,SVEC 
6 RIGOLS(I)»0.0 
C 

500 KXT*Hrm 

WRITE (8W. 1000) HIT 

C 

CALL SOLVES (AK ,V 1U ^BUBBKl) .WRCUMUmD*! ) .J3X1G . ISL.HDOF^HSOB) 
CALL ORTHO (VL.TR, HDOF, HRMOD, HVEC) 

C 

XJ*0 

DO 10 J*HRK0D«-1,HVEC 
C 

C CALCULATE THE UPPER PART OP AXX (SYMMETRIC) 

C 

DO 10 IaNRMODal.J 
TR*0.0 

DO 11 K-l.HDOF 

11 TR»TR0VL(K,I)*VR(K.J) 

IJ=IJH 
AKK(IJ)aTR 
10 CONTIHUE 
C 

CALL MULT (AM.VRd.HRMOD+D.VLCl.HRMOD+D.ISL.IDIAO.NSUB.NDpF) 

C 

IJaO 

DO 20 JaHRMOD+1 ,NVEC 
C 

C CALCULATE THE UPPER PART OF AMM (SYMMETRIC) 

C 

DO 20 IaHRMOD+1 . J 
TRaO.O 

DO 21 K»1,HD0F 

21 TR“TR+VL(K,I)*VR(X, J) 

, IJaIJ+1 
AMM(IJ)*TR 
20 CQHTIHUE 
C 

CALL JACOBI (AKX,AHH,XX t EIGY(HRMOD+l),NSMAX,TQLJAC,SSUB,NV) 

C 

C ORDER EIGENVALUES k EIGEHVECTORS 
C 

30 IS-0 

DO 40 I»HRH0D*1,HYEC1 

IF (EIGY(I*1) .GE.EIGV(I)) GO TO 40 
IS«1 

TR«EIGV(I*1) 

EIGV(I*1)=EIGY(I) 
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EIGV(I)=TR 
DO 41 J-i.HSUB 

TBaXX(J*(I"HRMOD)*HSUB) 

XX(J*(I-HRMOD)*HSUB)aXI(J>(I-HMOD-l)*HSTO) 

41 XX(J*(I-HRM0D-1)*HSOT)»TS 

40 CONTINUE 

IF (IS.EQ.l) GO TO 30 
C 

C SUBSPACE CONVERGENCE TEST 
C 

IC0HV»0 

DO 60 I»HRHOD'H,HVEC 

TRaDABS((EIGOLD(I)-EIGV(I) )/EIGV(I) ) 

El GOLD ( I ) =EIGV (I) 

EIGV(I)*TR 

IF (TR.GT.TOLEIG.AHD.I.LE.NEIG) IC0HV»1 
SO COHTIHUE 

WRITE (HW.1001) (EIGV(I) , I»1,NVEC) 

if (laarv.Eq.o) go to ioo 

c 

IF (HIT . LE . NITNAZ) GO TO 70 
WRITE (NV.1002) 

GO TO 100 
C 

C UPDATE EIGEN VECTORS 
C 

70 DO 80 Ial.NDQF 

DO 80 J-l.HSUB 
TR*0 . 0 • 

DO 81 K=1,HSUB 

81 TR=TR+VR(I,K*NRH0D)*XX(K*(J-1)*HSUB) 

VL ( I , J ♦NRMOD ) =TR 
80 COHTIHUE 

DO 90 I*1,HD0F 

DO 90 J*NRM0D+1 ,HVEC 
90 VR(I,J)»VL(I,J) 

GO TO 600 
C 

C CALCULATE FIHAL EIGENVECTORS 
C 

100 DO 110 I*1.HD0F 

DO 110 J»1,HSUB 
TR*0.0 

DO 111 K*1,NSUB 

111 TR*TR* VL(I , K+NRMOD ) * XX ( K* ( J- 1 ) »NSUB ) 

VR(I,J+HRHOD)»TR 

110 COHTIHUE 

DO 112 I»1,NVEC 

112 EIGV(I)aEIGOLD(I) 

C 

RETURN 

C 

1000 FORHAT ( SX , 12HITERATI0H HO, 15) 

1001 FORMAT (6(2X,1PE10.3)) 

1002 FORMAT (5X.24HWE ACCEPT CURRENT VALUES) 

1003 FORHAT (//20X, 'SUBSPACE ITERATION ROUTINE'//' NB OF EIGENVALUES*' , 
1 IS/' SB OF VECTOR* ',16/* HB OF DOF* ',15/' TOLERANCE* ' , 1PE10 . 3/) 

1004 FORMAT (• HB OF RIGID MODES** .IS//) 

C 
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END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE INVECT (AX, AM.VL, VR.IDIAG .NDOF.NVEC) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION AI(i),,AM(l),VL<l),VR(NDOF,NVIC) e imG(i) 

C 

hd-ndof/hvec 

c 

DO 10 :>i,NDQF 

neisusd) 

ra(i,i)«AM(xx) 

n.(i)»AM(ii)/«(ii) 

DO 10 J-2.NVEC 
VR(I,J)a0.O 
10 ccnmsuE 
c 

LLaNDOF-ND 

C 

DO 20 4*2»JREe 
TR=0 . 0 
C 

DO 30 I»1,LL 

IF (VL(I).LT.TR) GO TO 30 

TRsVUX) 

IJ«I 

30 CONTINUE 

C 

DO 40 I»LL,NDOF 

IF (VUl) .LE.TB) GO TO 40 
TRaVL(I) 

IJ»I 

40 CONTINUE 

C 

VL(IJ)»0.0 

LL»LL-ND 

VR(IJ,J)»1.0 

C • v 

20 CONTINUE 
C 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cc 
SUBROUTINE MULT (AM,VR,VL,ISL,IDIAG,HVEC,NDOF) 

IMPLICIT REAL*8 (A-H.Q-Z) 

DIMENSION AM(l),VR(NDOF,HVEC),VL(NDOF,HVEC),ISL(l) f IDIAG(l) 

DO 500 IVal.HVEC 
DO 100 I»1,ND0F 
TR=>0.0 
IJalDLAG(I) 

IK»ABS(ISL(I>) 

KK»I 

DO 110 X>IX,I 

TR=TR+AM(IJ)*VL(KK,IV) 

IJaIJ-1 

KK=KK-1 

110 CONTINUE 

C 

IF (I.Eq.NDOF) GO TO 99 
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IX«H-1 

DO 120 K^IK.HDOF 

I? (I.LT.ABS(ISL(K))) CO TO 120 
IJ«IDIiG(K)-K4.I 
TR*TR+AM(IJ)*VL(K , IV) 

120 CONTINUE 

99 VR(I , IV)*TR 

100 CO N T INUE 

500 COHTIHUE 

HETURH 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE JACOBI (AK,AH,XX,EIGV,BSHAX,TOL,N,NV) 

IHPLICIT REAL *8 (A-H.O-Z) 

DIMENSION AK(l),AM(l),XX(N,N).EIGV(l) 

C 

C INITIALIZE 
C 

ZEROoO.O 
TW0«2.0 
DO 10 I»1,H 

II*I*(I-1)/2*I 

IF (AK(II).LE. ZERO. OR. AK(II).LE. ZERO) GO TO 900 
EIGV(I)»AJC(II)/AM(II) 

DO 20 J°1 ,N 
XX(I,J)=ZERO 
20 CONTINUE 

XX(I.I)»1.0 
10 CONTINUE 
C 

C SET COUNTER 

C 

NSWEEP»0 

NR=N-1 

SOO NStfEEP*»NSVEEP+l 
C 

C CHECX IF ZEROOING IS REQUIRED 

C 

EPS=0 . 0 1**NSWEEP 
EPS=EPS*EPS 
DO 150 J»1,NR 
IIK»J*1 

DO 150 XallX.N 
JJ»J*(J-1)/2*J 
KK«K* (K- 1 ) /2*K 
JK=K*(K-l)/2+J 

EPSAK«( AK ( JK) * AK ( JK ) ) / ( AX ( J J ) * AK (KK) ) 
EPSAM«(AH(JK)*AM(JK))/(AH(JJ)*AM(KK)) 

IF (EPSAX.LT.EFS.AND.EPSAM.LT.EPS) GO TO 150 
C 

C CALCULATE ROTATION ELEMENTS 

C 

AKK»AK(KK)*AM( JK)-AM(KK)*AK(JK) 

AJJ»AK(JJ)*AX(JK)-AH(JJ)*AK(JK) 

AB=(AK(JJ)*AM(XX)-AX(XX)*AM(JJ))/TVO 

CHECK=AB*AB+AKX*AJJ 

IF (CHECK.LT. ZERO) GO TO 900 

SqCH=DSQRT( CHECK) 
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D1*AB+SQCH 

D2alB“SqCH 

DEH»D1 

IF (DABS(D2) .GT.DABS(Dl)) DEH=D2 
I? (S19.lS.ZZaO) GO TO 46 
CA«ZEBO 

CGs»AS(JI)/AX(KK) 

GO TO SO 

46 ClaAKX/DEH 

GGs-AJJ/DEH 


PMFtffil 6EHERALIZED BOTH I OH 


SO 


60 

70 


80 

90 


100 

110 


JPlaJ+t 

JM1»J-1 

KP1»K*1 

KH1»X~1 

IF (Jli.LT.l) GO TO 70 

00 60 

IJ s J*JMl/2+I 

IKaK*KM/2*I 

AKJaiK(IJ) 

AKX-AX(IK) 

AMJaAM(IJ) 

AXK»AM(XX) 

AX(IJ)=1XJ*CG*AXX 

AH(IJ)=AJU+CG*AHX 

AX(IK)=AXX*CA*AXJ 

AM(IK)»AHK«-CA*AJIJ 

C08TIBUE 

IF (KP1.GT.H) GO TO 90 
DO 80 I»XP1,H 
JI*»I*(I“l)/2+J 
KI»I*<I-1)/2*X 
AKJ«AK(JI) 

AMJaAM(JI) 

AKK 3 AK(KI) 

AKX 3 AH(XI) 

AK(JI) a AXJ+CG*AXX 

AH(JI) 3 AHJ*CG*AMK 

AX(XI) 3 AXX+CA*AXJ 

AM(XI) 3 AHK+CA*AHJ 

COHTXNUE 

IF (JP1.GT.XH1) GO TO 110 
DO 100 I-JP1.XH1 
JI«I*(I-1)/2*J 
IX*8>(X“1)/2*I 
AXJ»AX( JI) 

AMJ»AM(JI) 

AKK=AX(IX) 

AMKbAN(IK) 
AX(JI)«AXJ*CG*AXX 
AH( JI)»AXJ+CG*AKX 
AX(IK) 3 AKK+CA*AKJ 
AX(IK)*AHK+CA*AMJ 
CQBTIHUE 
AKX«AX(XX) 

IKK 3 AX (XX) 

AXJ a AX(JJ) 
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o o o o o o n ct o non 


AXJaAM(JJ) 

AK(KX)»1KKVTV0*CA'*AX(JK)+CA*CA*AKJ 

AM(KX)«AXX+>TVO*CA*AH(JX)*CA*CA*AXJ 

AX ( J J ) a AX J ♦TVO*CG* AX ( JX ) •*-CG*CG« AKK 

AM(JJ)»AMJVnfO*CG*AM(JK)*CG*CG*AKK 

AX (JX) “ZERO 

AN(JK)»ZERO 

UPDATE EIGEHVECTOR FOR THIS ROTATION 


120 

150 


DO 120 I»1,N 
XXJ«XX(I,J) 
XXK a XX(l,K) 
XX(I t J)aXXJ*CG*XXX 
XX(I f K) a XXK*CA*XXJ 
CONTINUE 
CONTINUE 


UPDATE ggEHVALDES * CHECK CONVERGENCE 

NT«N*(N-H)/2 
ICOHVaO 
DO 160 I»1,H 

II*I*(I-l)/2+I 

IF (AK(II).LE. ZERO. OR. AH(II).LE. ZERO) GO TO 900 
TR=AX(II)/AM(II) 

DEH» (TR-EIGV ( I) ) /TR 
EIGV(I)*TR 

IF (DABS(DEN) .GT.TOL) IC0NV=1 
160 CONTINUE 

IF (ICOHV.Eq.l) GO TO 499 

CHECK OFF DIAGONAL TERNS 

EPS*TQL*TOL 
DO 170 J»1,NR 
IIKaJ+1 

DO 170 K»IIK.H 
JJaJ*(J-l)/2+J 
KX*K*(K-l)/2+K 
JK»K*(K-l)/2+J 

EPSAK a (AK(JK)*AK(JK))/(AK(JJ)*AK(KK) ) 
EPSAN»<AH(JK)*AH(JK))/(AM(JJ)«AX(KK)) 

IF (EPSAK.LT.EPS.AND.EPSAN.LT.EPS) GO TO 170 
GO TO 499 
170 CONTINUE 

SCALE EIGENVECTORS 


179 DO 180 I»1,N 

II»I*(I-l)/2+I 

AXK*DSQRT(AM(II)) 

DO 180 Jal.N 

XX(J,I)»XX(J,I)/AKK 

180 CONTINUE 

* 

RETURN 

* 

499 IF (NSVEEP.LE.NSNAX) GO TO 500 
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WRITE (NV.1000) 

GO TO 179 
C 

900 WRITE (IV, 1001) 

STOP 

C 

1000 FORMAT (6X,34BM0 CONVERGENCE AT BSMAX ITERATIONS) 

1001 FORMAT (SX.46HERR0R IB JACOBI : XATRII HOT POSITIVE DEFINITE) 



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SOiaOOTIBE SOLVES (AI.VL.VR.IDIAfi.ISL.BDOF.NVEC) 

C ? 

C PURPOSE : SOLVES THE SINGULAR PROBLEM : AX x VL ■ VR 
C TEE SINGULAR COLUMNS INTO AX ARE INDEXED BY 

C TEE NEGATIVE VALUES OF ISL, TEE CORRESPONDING 

C TERMS OF THE SOLUTION ARE PUT TO 0 

C 

c 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION AX(l) ,VL(NDOF,NVEC) ,VR(NDOF,NVEC) ,IDIAG(l) ,ISL(1) 

C 

DO SOO IV»1,HVEC 

DO 60 I»1,ND0F 
60 VL(I,IV)»VR(I,IV) 

BACKSUBSTITUTE 

DO 100 XC>2,HD0F 
TRsO.O 
ICl-IC-1 

IH1«IDIAG(IC)-IC 
IKalSL(IC) 

IF (IX.LE.O) THEN 
VL(IC,IV)aO. 

GOTO 100 
END IF 

IF (IK.GT.IC1) GO TO 100 
DO 120 X-IX.IC1 

TR*TR+AK( IN1+K)*VL(K , IV) 

CONTINUE 

VL(IC,IV)=>VL(IC,IV)-TR 
CONTINUE 

SOLVE DU»U 

DO ISO ICal.HDOF 

VL(IC.IV)3VL(IC,IV)/AK(IDIAG(IC)) 

ISO CONTINUE 

BAK SUBSTITUTE 

IICaBDOF 
DO 200 IC»2,HD0F 
TRaVL(IIC.IV) 

IC1»IIC-1 
IXalSL(IIC) 
nri»iDiAG(iic)-iic 
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IF ((IK.GT.IC1).0R.(IX.LE.05) GO TO 221 
DO 220 K*IK,IC1 

VL(K ,IV)=VL(K,rV)-lX(ISl+K)*TR * 

220 CONTINUE 

221 IICaIIC-1 
200 COBTXSUB 
BOO CONTINUE 

RETURN 

EHD 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE FACT (AX,IDIAG.ISL,NDOF,HV) 

c 

C PURPOSE : LDLT DECOMPOSITION OP THE POSITIVE SERI-DEFINITE 
C MATRIX AX, THE SINGULAR COLUMNS OF AX ARE INDEXED 

C BY A NEGATIVE VALUE OF IDIAG> 

C 

C 

UEELICTT B£AL«8 (A-H.0-Z3 
DIMENSION AX(1),2DIAG(1),ISL(1) 

C 

C DETERMINE MIN A MAX 
C 

trsdabs(akcD) 

AMIN-TR 

AMAX*TR 

DO 10 IL»1,ND0F 

TR»DABS(AK(IDIAG(IL))) 

IF (TR.LT.AMIN) AMIN“TR 
IF (TR.GT.AMAZ) AHAX»TR 
10 CONTINUE 

ZER0 3 (AMAZ+'AMIH)*1 .0D-10 

c 

C LOOP OVER COLUMN 
C 

DO 100 ICal.NDOF 
MIC=ISL(IC) 

IClaIC-1 

MIC1=HIC-*T 

IN2*IDIAG(IC)-IC 

IF ((HIC.LT. 1) .OR. (MIC.GT. IC) ) GO TO 901 
C 

C CALCULATE GS 

C 

IF (MIC1.GT.IC1) GO TO 150 
DO 120 ILoHICl'ICl 

IF (XDIAG(IL) .LT.O) THEN 
AX(IN2+IL) 3 0 
GOTO 120 
END IF 

MILalSL(IL) 

ILl=IL-i 

MIH»MAXO(MIL.MIC) 

IH1*IDIAG(IL)-IL 

TRaO.O 

IF (MIM.GT.IL1) GO TO 120 
DO 130 KaMIM.ILl 

TR=TR+AK(I51+K)»AXCIS2+K) 

130 CONTINUE 
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IN»XH2*IL 
AK(XM)*AK(IN)-TR 
120 CONTINUE 

C 

C CALCULATE LAO 

e 

ISO TR»0.0 

IF (IIC.6T.IC1) 60 TO 201 
DO 200 ILaMIC.ICl 

IF (IDIAG(IL).LT.O) GOTO 200 
AG*AX(XS2*IL) 

AL*AG/AX(XDIAG(IL)) 

AK(IH2*IL)«AL 

TR“TR+AL*AG 

200 continue 

201 IHsIDIAaClC) 

AK(XN)*AX(IN)-TR 

IF (AK(IN) .LT.ZERO) IDIAG<IC)— IDXAG(IC) 

ioo co st un i s 
RETURN 

901 WRITE (NV.lOOl) 

STOP 

1001 FORMAT (SX,29H«**ST0P ERROR IN ID I AG VECTOR) 

1010 FORMAT (SI, ’ CONDITIONING OF THE STIFFNESS MATRIX '/2X, 

1 'KIN D1AG TERMs' .1PE10.3,' XAZ DIAG TERH=' .E10.3) 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE NULL (AK,AH s VL,VR,IDIAG,ISL,NDOF,NRHOD) 

C 

C PURPOSE : CALCULATE THE NULL SPACE OF AX AND PUT AN 
C ORTHONORMALXSES BASE OF THIS SPACE INTO THE 

C NRMQD FIRST VECTORS OF VR. 

C VL ■ AM x VR AFTER EXECUTION 

C 

C 

IMPLICIT REAL*8 (A-H.Q-Z) 

DIMENSION AK(1) ,AM(1) ,IDIAG(1) ,ISL(1) ,VL(NDOF,*) ,VR(NDOF,*) 

C 

C STORE THE SINGULAR COLUMNS INTO THE BEGINNING OF VR 
C 

C THE SINGULAR EQUATION ARE NOV INDEXED BY NEGATIVE VALUES 
C INTO ISL INSTEAD OF ID I AG 
C 

NRMQD >0 

DO 1 IC*1,ND0F 

IF (IDIAG(IC).GT.O) GOTO 1 
IDIAG(IC)a-IDXAG(XC) 

NRHOD*NRMQD+ 1 
MIC*ISL(IC) 

IH=IDIAG(IC)-IC 
DO 2 K»i,MIC-l 

2 VR(K,NRM0D)=0. 

DO 3 K-MIC.IC-l 

VR(X,NRMOD)»AX(IN+K) 

AK(IN*K)»0. 

3 CONTINUE 
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ISL(IC)«-ISL(IC) 

Va(IC,SHBOD)»-l. 

AK(XDXAG(IC)) a l. 

DO 4 X«ICn,HDOF 
4 VR(K ,HRHQD)*0 . 

1 COHTXHUE 
C 

C BAXSUBSTITUTE 

e 

DO 200 Hal.HBXOD 
IIC-HDOF 

DO 200 XC*2,HD0F 
TRaVR(IIC.H) 

IClaIIC-1 

IE«ISL(IIC) 

m»IDIAO(IIC)-IIC 

IF ((IX.GT.IC1) . OE. (IK.LE.O)) GO TO 221 
DO 220 K-XK.ICl 

220 COHTXHOB 

221 IICaIIC-1 
200 COHTZHOE 

C 

C ORTHOGOHALISATIOH , ; ; 

C 

DO .10 H»1,HRH0D 
DO 20 Kal.H-l 
TRaO . 

DO 30 lal.HDOF 

30 TR»TR*VI.(I,X)*VRU 1 H) 

DO 40 lal.HDOF , 

40 VR(I ,H)aVR(I.,H)“TRaVR(I ,X) 

20 COHTIHUE 

CALL HULT(AH,VL(1 ,H) ,VR(1 ,H) .ISL.IDIAG , 1 ,NDOF) 

TRaO. 

DO 50 lal.HDOF 

SO TRaTRm<I,H)*VL(I-,H) 

TRal/SQRT(TR) 

DO 60 Xal.HDOF 

VR(I,H)aVR(I.H)*TR 
VL(I,H)»VL(I,H)*TR 
60 COHTIHUE 

10 COHTIHUE 

RETURH 

EHD 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTXHE ORTHO ( VL, VR, HDOF, NRMOD ,HVEC) 

C 

C PURPOSE : ORTHOGOHALXSE THE HSOB LAST COLUKHS OF VL (LAST 
C EVALUATED SOLOTIOH) WITH RESPECT TO THE HULL SPACE 

C OF A. FOR THE AH SCALAR PRODUCT 

C 
C 

IMPLICIT REAL*8 (A-H.O-Z) , 

XHTEGER HDOF.HRHOD.HVEC 

DIKEHSIOH VLCNDOF ,HVEC) , VR(HDOF ,HVEC) 

HSUBaHVEC-HRMOD 
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DO 1 Jal.HSUB 
DO 2 2*i,HRSO0 
S-O. 

DO 3 Kat.IDOF 

3 S»S*VL(K , I) *VL(K .NRHOD* J ) 

DO 4 L»l t HDOF 

4 7L(L,HHJI0D+J)aVI,(L,HIUC0D+J)-S*Va(L,I) 

a coHTuroE 

1 COSTXHUE 

KSTOBV 

END 


ccccccccccecccccccccceccccccccccccccccccccccceccccccccccccccccccccccc 


File: animout.f 


C-Module ANXMOUT 

OAutbor K. Alvin 

C=Date June 1990 

C^Bloci; Fortran 

C 

C subroutine AHIHOOT 

C 

C Purpose: 

C This subroutine produces an output file to be used to 

C visualize tbe simulation using MESH. 

C 

C . 

c 

subroutine AHMOOT(q I id,nnp .time ,out) 
real*8 q(l) .time ,v(3) 
integer out,id(6,<),uup,i,k 
vrite(out,'(gl5.S)’) time 
do 100 isl.nnp 
do 200 k«l,3 

if (id(k,i) .ne. 0) then 
v(k) » q(id(k,i)) 
else 

v(k) » 0. 
end if 

200 continue 

write (out, 1000) (v(k),k*l,3) 

100 continue 

1000 format(3glS.S) 
return 
end 


File: stiffrc.f 


C subroutine STIFFRC ( r .fact ,kq) 
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racoraiva aubroutiaa STIFFRC(v,f act ,kq) 
rad*8 T(l),fact,kq(l) 
includa 'sharad.iac' 

C ASSEMBLE EACH ELEMEHT MASS AMD STIFFNESS 

do 100 ii*l .adoaain 

CVOt CHCALL 

do 110 jjal.aald(ii) 
a a alau»(jj,ii) 
call ELEFRC(7,*act,kq,a) 

110 eoatiaaa 

100 eoatiaaa 

ratora 

•ad 

C aabroatiao ELEFRC(v,tact ,kq,n) 

racuraiva anbrontiaa ELEFRC(7,f act ,kq,a) 

raal*8 7(1) ,fact,kq(l) 
iatagar a 

iaduda ’ aharad.iac' 

C LOCAL VARIABLES 

para»atar(HAXSEQ°24) 
raal*8 ak(MAISEQ .MAISEQ) 
iatagar lm(HAXSEQ) t aaaq 

do 20 k»l,4 
jaiz(k.a) 

it ((atypa(a) .aq.l) ,aad.(k.gt .2)) j » 0 
do 30 ial,6 
kkafla(k-l) ♦ i 
it (j .aa. 0) thaa 
la(kk) » idCi.j) 
alaa 

la(kk) > 0 
•adit 

30 eoatiaaa 

20 eoatiaaa 

acaqsl2 

call L0A0SK(ak,a,aaaq) 

call ESTIFVMCak.la.asaq.T.kq.fact) 

ratura 

•ad 
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C subroutine LOADSK(ak,a,naeq) 

recursive subroutine LOADSKCsk.n.nseq) 

include ’ shared. inc’ 

6 

real*8 sk(nseq,l) 
integer a, nseq 

k«0 

do 10 jal t nseq 
do 20 i«i,j 

k«kn 

ek(i, j)aestifm(k,n) 
sk(j.i)»sk(i.j) 

20 continue 

10 continue 

return 

end 


File: estifvm.f 


C subroutine ESTIFVH(ak,lni,aaeq, v.kq.fact) 

recursive subroutine ESTIFVH ( ak , la , na aq , v , kq , f act ) 
C ARGOKEHYS 

real«8 sk(nseq.l) ,v(l) ,kq(l) ,iact 
integer ln(l),nseq 

do 20 j • 1, nseq 
k ■ l*(j) 

if (k . eq. 0) goto 20 
do 10 i = 1, nseq 
m a lffi(i) 

if (a .eq. 0) goto 10 
kq(m) a kq(m) ♦ sk(i, j)*v(k)*fact 
10 continue 

20 continue 

return 

end 

CaEnd Fortran 


File: renum.f 


CaDECK REHUH 

C=PURP0SE- REHUHBERS THE GRID P0ZHTS TO KI51HIZE PROFILE STORAGE 
C=A0TH0a V S BELYZH and DOC HGU7EH 7-5-90 
C 

subroutine REHTJH 
C 

include ' shared. inc’ 
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Initial!*# vector 
maxtry*2*nap/ 3 ♦! 


ater»s»nnp*aaxtry 
do 22 j*i,nteras 
22 iadjcy(j)«0 

do 10 i*l,nap 
10 icoaat(i)«0 

do 1 i*l,nel 
nodea»ix(l,i) 
aodeb>ix(2,i) 

if ((aodoa .eq. 0).or.(aodeb . sq. 0)) go to 1 
ieount (nodea) «icouat (aodoa) *1 
ieount (nodab) *icouat (nodab ) > 1 
lamicouatfaodoa) 
ib*icount (nodab) 

if (ia.gt .maxtry .or. ib.gt. maxtry) go to 34S 
locate a (aodea-l)*maxtry+ia 
ladjcy(locate)*nodeb 
loeatea(nodab-l)*maxtry*ib 
iad jcy (locate) *nodaa 
1 contiaua 


ii a 0 

do 37 i*l,ntarma 

if ( iadjcy(i) .eq. 0 ) go to 37 
ii«ii+l 

iadjcyCii)aiadjcy(i) 

37 coatiaaa 


last*0 

do 2 i*l,nnp 

laat =las t ♦icoont (i ) 

2 coatiaaa 

jjoicooat(l) 
icouat(l)«l 
do 3 i*2 ,nnp*l 
kk«icount(i) 
i cooat ( i ) aicooat ( i- 1 ) ♦ j j 
jj»kk 

3 coatiaaa 

go to 666 
346 vrite(6,&55) 

566 format(2x, • arror in dimansion for KAZTE7 !! ') 
556 coatiaaa 
c ««*****«******«**a«*aaa 

call GZKRClKnnp.icoant.iadjcy.para.mask.xls) 

ratoxa 

and 

cxxmxxxxxxxxxmxx 


subroutine ganrca(naqns , xad j , ad j acy , porm .mask , zls ) 

c rafaraaca: computer solution of large sparse positive definite 

c systems, alaa gaorga ft josepb w-h lin 

c (prentive-hall.iac. .eaglesood cliffs, HI 07632) 

integer adjncy(l) ,aask(l) ,parm(l) ,zls(l) 
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iateger xadj (1) , ccsize, i.neqna ,nlvl,aus,root 
do 100 ial,neqns 

100 coatiaao 
nasal 

do 200 1*1 ,neqne 

i£( aaak(i).eq.0 ) go to 200 
rootsi 

call fnroot(root , xadj , ad jncy, mask, nlvl.xls ,perm(num) ) 
call rca(root, xadj .adjacy .mask, perm(aum) , ccaize.xls) 
aasenuBtccsize 
if(au».gt .neqae) go to 987 


200 coatiaao 

c pera(nev node) sold aode 

c now, aaakCold node)* nee aode 


987 coatiaao 

do 11 aew»l, neqae 
iold*pexs(aev) 

®&sk(iold)*ne* 

e arite(6,*) ’ iold,mask(iold) * ' ,iold,mask(iold) 

11 coatiaao 
retara 
ead 

cxxxxx.,zm%xxxxmxxxxxmxxxxmxxxxxxxxxxxxxxxxxxxxx 

subroutine faroot (root , xadj .adjncy, mask, nlvl.xls, Is) 
integer adjacy(l) ,la(l) ,mask(l) ,xls(l) 
iateger xadj(l) , ccsize, j , jatrt,k,kstop,kstrt, 

$ misdeg.nabor.ndeg.nlvl, node, nualvl, root 

call rootla(root,xadj , ad j acy , mask , alvl , xla ,1s) 
ceaizesxla(alvl+l)-»l 

if(nlvl.eq.l .or. alvl.eq. ccsize) return 
100 jatrt»xla(nlvl) 
miadegaccaixe 
root»la(j8trt) 

if ( ccsize. eq, jatrt) go to 400 
do 300 jsjstrt, ccsize 
nod*«Xs(j) 
adegaO 

kstrt»xad j (node) 
katopaxadj (node+l)-l 
do 200 kakatrt .kstop 
aabor=adjncy(k) 

if( mask(aabor) .gt. 0) ndeg=ndeg+l 
200 coatiaao 

if (adeg.ge.aindeg) go to 300 
root*node 
nindeg*ndeg 
300 coatiaao 

400 call rootla (root, xadj , ad j acy, mask, nualvl.xls ,1s) 
if (aaalvl.le.alvl) return 
alvl»nanlvl 

if (alvl.lt. ccaize) go to 100 

retara 

ead 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

aabroatiae rca(root,xadj .adjncy, mask, perm, ecaize.deg) 
iateger adjncy(l) ,deg(l) ,aask(l) ,perm(l) 
integer xadj(l) ,ccaize,fabr,i,j ,jstop,jstrt,k,l,lbegin, 
8 lnbr.lperm.lvlead.abr, node, root 
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call degr**(root,x&dj .adjney, mask, d*g,ecaixe, pens) 

aa*k(root)«0 

if (ccsixe.le.l) return 

lTl*nd*0 

lnbr*l 

100 lbagiaalvlandal 

lTl«ad«bbr 

do 600 i<*lb*gin,lvl*ad 
node*perm(i) 
j strt*xad j (nod*) 
j»top*xadj (nod*+l)-l 
fnbr*labr*l 
do 200 j*jstrt , jatop 
abr»adjney(j) 

if (maak(nbr) . eq.0) go to 200 

labr«labr*l 

Bask(nbr)>0 

p«ra( lnbr) ®nbr 

200 continue 

if(fnbr .g*.lnbr) go to 600 
k»£nbr 

300 l»k 

k»k*l 

nbr»pora(k) 

400 if(l.lt.fnbr) go to S00 

lp*ra*p*rm(l) 

if( d«g(lp«rm) . 1* . dog(nbr) ) go to S00 
p*rm(l*l)*lp*ra 

lBl_i , 
go to 400 

SOO p*na(l*l)°nbr 

if (k. It. lnbr) go to 300 

600 coot inn* 

if (lnbr.gt .lvlend) go to 100 
k*ccsiz*/2 
l*ccsixe 
do 700 i*l,k 
lp*rw*p*ra(l) 

P*na( 1) *p*r»( i) 

p*xm(i)*lp*ra 

1 » 1-1 

700 continue 

r*tnm 
•nd 

subroutine rootls(root ,xsdj , ad jncy.aaak, nlvl.xls ,1s) 
int*g*r adjncy(l) ,ls(l) ,mask(l) ,xls(l) 
integer xadj(l) ,i,j , jstop, jetrt.lbegia.ccsize.lvlend, 
$ lxsine.nbr ,alvl, node, root 

mask ( root )<*0 
la(l)*root 

nlrl»0 

lxl*nd«0 

eesiXMl 

200 lbegin*lvl*nd>l 
lrl«nd>eceize 
nlvlenlul+i 
xls ( nlvl) »lb*gin 
do 400 i*lb*gin, lvlend 
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aedasls(i) 
j strtax&dj (node) 
j stop*xad j (aod«> i ) -i 
12 (jatop.lt. jatrt) go to 400 
do 300 j* jatrt, jatop 
abroad jncy(j) 

if ( aaak(nbr) .aq.O) go to 300 
eeaizaseeaiza+l 
la(ceaiza)aabr 
B&sk(nbr)*0 
300 eoatiaaa 
400 eoatiaaa 

lvsiz**ccsiza-lvlaad 
if (lvaiza. gt .0) go to 200 
xls(nlvl*l)alvlead+l 
do 800 i»l,ecsiza 
noda«la(i) 
aaak(noda)*l 
500 Croatian* 

ratara 
and 

c mxzmmmama 

subroutisa dagr aa (root, xad j .adjncy.naak.deg.ccsiza.ls) 
intagar adjncy(l) ,dag(l) ,la(l) .oask(l) 
iatagar xadj(l) ,ccaiza,i,idag, j , jatop, jstrt, 

$ Ibagin , lvland , lvaiza . abr , node . root 

la(l)«root 

xad j (root ) "-xad j (root ) 

lvlaadaO 

ccaiza»l 

100 lbagia»lvland>l 
lvlaadaceaiza 
do 400 islbagia, lvland 
aoda»la(i) 
j a t rt =-xad j ( noda ) 
j»top*iab»( xadj(noda+l) ) -1 
idagaO 

if (jatop. It. jatrt) go to 300 
do 200 j® jatrt, jatop 
nbr®adjney(j) 

if( naak(abr) .aq.O ) go to 200 
idag*idag+l 

if (xad j (abr) .It. 0) go to 200 
xad j (nbr)a-xadj (abr) 
ccaiza*ccaiza+l 
la(ccaiza)®nbr 
200 eoatiaaa 

300 dag(noda)=idag 
400 eoatiaaa " 

lvaiza*ceaiza- lvland 
if (lvaiza. gt.0) go to 100 
do 500 i*l,ccaize 
noda®la(i) 

xadj (noda)s-zadj (node) 

500 eoatiaaa 

retarn 
and 


File: kfilter.f 
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C subroutine KFILTER 

recursive subroutine KFILTER 
include ’shared.inc* 

f 

do 10 i » l.ndof 

go(i) * delsq*(f (i)+go(i))+delta*pe(i)4aass(jdiag(i))*qe(i) 
gk(i) « delta*(f (i)*gk(i))+pe(i) 

10 continue 

call SOLYER(eo,go,jdiag,ndof ,2) 

Activate EBE computations for internal force by using STIFF8.C 
subroutine. Otherwise use PKVMAD (profile matriz/vector mult-add 

C call ?ff¥KAD(stii .jdiag, go, ndaf, -delta, gk,A.dO) 

call STIFFRC (go, -delta ,gk) 

do 100 i a l.ndof 

qe(i) o 2.*go(i) - qe(i) ' 

pe(i) » 2.*gk(i) - pe(i) 

100 continue 

return 

end 


<3 
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